{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thermal Transport from NEMD\n",
    "# 1. Introduction\n",
    "- In this tutorial, we use the NEMD method to study heat transport in graphene at 300 K and zero pressure. Here we consider the ballistic regime. The diffusive regime can be found in the HNEMD example. The spectral decomposition method as described in [[Fan 2019]](https://doi.org/10.1103/PhysRevB.99.064308) is used here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing Relevant Functions\n",
    "- The inputs/outputs for GPUMD are processed using the [Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/) and the [thermo](https://github.com/AlexGabourie/thermo) package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pylab import *\n",
    "from thermo.gpumd.io import ase_atoms_to_gpumd\n",
    "from thermo.gpumd.preproc import add_group_by_position\n",
    "from thermo.gpumd.data import load_shc,load_compute\n",
    "from ase.build import graphene_nanoribbon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Preparing the Inputs\n",
    "- We consider a graphene sheet of size of about 25 nm x 43 nm (40400 atoms). The transport is in the $y$ direction. We divide the length in the $y$ direction into 9 groups. Group 0 (a small one) will be fixed. Group 1 (about 9 nm long) will act as a heat source and group 8 (about 9 nm long) will act as a sink. The remaining middle part is evenly divided into a few groups (group 2 to group 7).\n",
    "- We use a Tersoff potential [[Tersoff 1989]](https://doi.org/10.1103/PhysRevB.39.5566) parameterized by Lindsay *et al.* [[Lindsay 2010]](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.81.205441).\n",
    "\n",
    "## Generate the  [xyz.in](https://gpumd.zheyongfan.org/index.php/The_xyz.in_input_file) file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Atoms(symbols='C40400', pbc=[True, True, False], cell=[245.95121467478057, 430.26, 3.35])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gnr = graphene_nanoribbon(100, 101, type='armchair', sheet=True, vacuum=3.35/2)\n",
    "# Rotate to match example\n",
    "gnr.euler_rotate(theta=90)\n",
    "l = gnr.cell.lengths()\n",
    "gnr.cell = gnr.cell.new((l[0], l[2], l[1]))\n",
    "l = l[2]\n",
    "gnr.center()\n",
    "gnr.pbc = [True, True, False]\n",
    "gnr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Add Groups for NEMD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 3.9026, 89.5546, 132.1806, 174.8066, 217.4326, 260.0586, 302.6846, 345.31059999999997, 430.26]\n"
     ]
    }
   ],
   "source": [
    "split = array([l/100] + [l/5] + [l/10]*6)-0.4\n",
    "split = [0] + list(cumsum(split)) + [l]\n",
    "print(split)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Atoms per group: [400, 8000, 4000, 4000, 4000, 4000, 4000, 4000, 8000]\n",
      "Total atoms: 40400\n"
     ]
    }
   ],
   "source": [
    "ncounts = add_group_by_position(split, gnr, direction='y')\n",
    "print(\"Atoms per group:\", ncounts)\n",
    "print(\"Total atoms:\", sum(ncounts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ase_atoms_to_gpumd(gnr, M=3, cutoff=2.1, sort_key='group', order=range(9), group_index=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The <code>run.in</code> file:\n",
    "The <code>run.in</code> input file is given below:<br>\n",
    "```\n",
    "potential    potentials/tersoff/Graphene_Lindsay_2010_modified.txt 0\n",
    "velocity     300\n",
    "\n",
    "ensemble     nvt_ber 300 300 0.01\n",
    "fix          0\n",
    "time_step    1 \n",
    "dump_thermo  1000        \n",
    "run          100000\n",
    "\n",
    "ensemble     heat_lan 300 100 10 1 8 \n",
    "fix          0\n",
    "compute      0 10 100 temperature\n",
    "compute_shc  2 250 1 1000 400.0 group 0 4\n",
    "run          1000000\n",
    "```\n",
    "\n",
    "- The first line uses the [potential](https://gpumd.zheyongfan.org/index.php/The_potential_keyword) keyword to define the potential to be used, which is specified in the file [Graphene_Lindsay_2010_modified.txt](https://github.com/brucefan1983/GPUMD/blob/master/potentials/tersoff/Graphene_Lindsay_2010_modified.txt).\n",
    "\n",
    "- The second line uses the [velocity](https://gpumd.zheyongfan.org/index.php/The_velocity_keyword) keyword and sets the velocities to be initialized with a temperature of 300 K. \n",
    "\n",
    "- There are two runs. The first [run](https://gpumd.zheyongfan.org/index.php/The_run_keyword) serves as the equilibration stage.\n",
    "  - Here, the NVT [ensemble](https://gpumd.zheyongfan.org/index.php/The_ensemble_keyword) (the Berendsen thermostat) is used. The target temperature is 300 K and the coupling constant is 0.01 (dimensionless) for the thermostat. \n",
    "  - The [time_step](https://gpumd.zheyongfan.org/index.php/The_time_step_keyword) for integration is 1 fs. \n",
    "  - Atoms in group 0 will be fixed. \n",
    "  - The thermodynamic quantities will be output every 1000 steps. \n",
    "  - There are $10^5$ steps (100 ps) for this run. \n",
    "- The second [run](https://gpumd.zheyongfan.org/index.php/The_run_keyword) is for production. \n",
    "  - Here, two local Langevin thermostats ([ensemble](https://gpumd.zheyongfan.org/index.php/The_ensemble_keyword)) are applied to generate a nonequilibrium heat current. The time parameter in the Langevin thermostats is 0.1 ps (100 time steps). The target temperature of the heat source is $300+10=310$ K and the target temperature of the heat sink is $300-10=290$ K. Atoms in group 1 are in the heat source region and atoms in group 8 are in the heat sink region.\n",
    "  - Atoms in group 0 will be fixed. \n",
    "  - The [compute](https://gpumd.zheyongfan.org/index.php/The_compute_keyword) is used to compute the group temperatures and energy transfer rate. Here, grouping method 0 is used, and the relevant data are sampled every 10 time steps and averaged for every 100 data points before written out. \n",
    "  - The line with the [compute_shc](https://gpumd.zheyongfan.org/index.php/The_compute_shc_keyword) keyword is used to compute the spectral heat current (SHC). The SHC will be calculated for group <code>4</code> in grouping method <code>0</code>. The relevant data will be sampled every 2 steps and the maximum correlation time is $250 \\times 2 \\times 1~{\\rm fs} = 500~{\\rm fs}$. The transport directions is <code>1</code> ($y$ direction). The number of frequency points is 1000 and the maximum angular frequncy is 400 THz.\n",
    "  - There are $10^6$ steps (1 ns) in the production [run](https://gpumd.zheyongfan.org/index.php/The_run_keyword). This is just an example. To get more accurate results, we suggest you use $10^7$ steps (10 ns)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Results and Discussion\n",
    "### Computation Time\n",
    "- Using a laptop with a GeForce RTX 2070 GPU, the NEMD simulation takes about 20 minutes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure Properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "aw = 2\n",
    "fs = 16\n",
    "font = {'size'   : fs}\n",
    "matplotlib.rc('font', **font)\n",
    "matplotlib.rc('axes' , linewidth=aw)\n",
    "\n",
    "def set_fig_properties(ax_list):\n",
    "    tl = 8\n",
    "    tw = 2\n",
    "    tlm = 4\n",
    "    \n",
    "    for ax in ax_list:\n",
    "        ax.tick_params(which='major', length=tl, width=tw)\n",
    "        ax.tick_params(which='minor', length=tlm, width=tw)\n",
    "        ax.tick_params(axis='both', direction='in', right=True, top=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot NEMD Results\n",
    " - The [compute.out](https://gpumd.zheyongfan.org/index.php/The_compute_keyword) output file is loaded and processed to create the following figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['Ein', 'Eout', 'm', 'T'])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute = load_compute(['T'])\n",
    "compute.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = compute['T']\n",
    "Ein = compute['Ein']\n",
    "Eout = compute['Eout']\n",
    "ndata = T.shape[0]\n",
    "temp_ave = mean(T[int(ndata/2)+1:, 1:], axis=0)\n",
    "\n",
    "dt = 0.001  # ps \n",
    "Ns = 1000  # Sample interval\n",
    "t = dt*np.arange(1,ndata+1) * Ns/1000  # ns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArsAAAFTCAYAAAA5nMTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3xUVfrH8c9JJwlJSGhJ6F2kCBJABRRw7SJYsa2uZd21rG5Xd3/b1HXdVdfV1d21rGsFK4q9UIOCVEF67zUJSUhIn/P7Y4aQSYGQzMydmXzfr5evkHvv3PsEk4cnZ855jrHWIiIiIiISjiKcDkBERERExF9U7IqIiIhI2FKxKyIiIiJhS8WuiIiIiIQtFbsiIiIiErZU7IqIiIhI2FKxKyIiIiJhS8WutAjGmEhjzGpjzLQmvj7JGJNnjPmrr2MTEQk39eVcY8xsY0yjm/sbY35rjDlkjGnvnyilpVCxKy3FLUA/4I9NebG1thD4O3CXMaaLLwMTEQlDzcq5Hk8CFcAffBGQtFwqdiXsGWOigP8DvrTWftuMW/0T98/MfT4JTEQkDPkq53oGGZ4HbjXGdPJVfNLyqNiVluAiIBN4rTk3sdYeBD4GrjXGJPoiMBGRMOSTnOvxGhAF3OyDe0kLpWJXWoIbARfgNV/XGHOqMeZpY8wqY0yhMabYGLPUGHO7McY0cK+3gdbA5f4NWUQkZN1IPTn3CGNMK2PMY8aYncaYUmPMcmPM9fVda61dDmz03FOkSVTsSlgzxkQAZwKrPW+J1XQrcAmwHPg38DKQDDyNe35ufeZ7Po7zfbQiIqHtODn3iLeAy4A3cU9TSAdeNsb8tIHr5wPdjDE9fB2vtAwqdiXcnQSkAEvqOfdnoIu19hpr7a+stT8G+gKf4V6I1rX2C6y1m4GDwBl+jFlEJFQdK+ce0R0YYK39mbX2TuAUYD/wsDEmo57rF3s+Ku9Kk6jYlXB3ZFHDvtonrLXbrbWuWscqgWdx/2yMbeCe+2rcV0REjmow59bwkLW26Mgn1trdwD+AWOCqeq4/ci/lXWkSFbsS7lI9H/NrnzDGxBpjfmmMWezp5Wg9PSDf8VyS3sA984AYY0xrP8QrIhLKGsy5NWTXc2ye5+Pges7leT62bWpQ0rJFOR2AiJ+Vej7G1XPuXeACYC3wOnAAqAS6ATfgHmWoTyvAAiW+DFREJAwcK+cesb+eY0dGb5PrOdfK8/FwU4OSlk3FroS7A56PqTUPGmOycBe6nwIX1pzOYIy5Cnex25BUIN8z5UFERI6qN+fW0h7YUetYB8/HgnquP3KvA/WcEzkuTWOQcLcK9yhs71rHe3o+flR73i7HWARhjEnAPW/sO59FKCISPhrKuTWNrufYKM/H5fWc6+v5qLwrTaJiV8KaZyOIlUBWrVPbPR+9CltjzEjgh8e45alAJDDHVzGKiISLY+Tcmn5Tc2MeY0w6cDdQhrsdWW0jgHJggQ9DlRZExa60BO8DqcaYoTWOfYO7nc1kY8wsY8xfjTFvA3OBD49xr7M9H9/zT6giIiGvvpxb0xZgpWdjiadwj+a2B+6z1u6qeaGnKB4JfGat1ToJaZKAFrvGmHONMTONMXuNMWWe3VPeNMb0r3FNJ2PMU8aY+caYw54V8t0auF+cMeZvxpg9xpgSz2vGBOrrkZDxPO7dfK47csBaW4V7S8uXcL9Fdifu3o83Av88xr2uAZZYa5f6K1gRJzQ29yrvSiPUybm1XIG7681VuN9J2wfcYK2tbzOfSbgXqD3rhzilhTDW2sA9zJirgaG4R9UOAF2Ae4HOwEBr7TZjzFnAG7gbUkcC5wDdrbVb67nfa8CFwC+BzcAdwPnAadbab/399UjoMMa8BYwBujV1dMAYcyYwG7jOWuuLPd9FgkZjc6/yrjSGL3Ku5z6zcLeBPNkzSCFywgJa7NYbgDF9cbd++oW19jFjTMSRBUPGmFuA56g/4Q4GvgVusta+6DkWhXty/Dpr7YQAfhkS5IwxfXB/b/zSWvtEE+8xA2gDnGqd/sER8bHG5F7lXWksH+Xc0binlk2y1mrqmDRZMMzZzfV8rASoZ2V8QyYAFbhHIvC8thKYCpxrjGmoR6q0QNba9binKJQe59J6GWOScCfdW1XoSjhqZO5V3pVGaW7O9WiDeyBMha40iyN9do0xkbjfJusK/AXYC0w5wducDGyx1tZuMr0KiAF6ef5c87kqUgRjzL+a8fI/GGN8FouEFmttS/+ff0J5VzlXoNk5F2PMo76KRUKLr3KuU5tKfIO7hRPARmCctba+HVWOJRU4WM/xvBrnRUTEd5R3RSTkOFXsXg8kAT2AXwBfGGNG1bcIzR9O/t2nTBySwa2je9A1LSEQjzyuI6OFwfYOebDGBYqtKYI1LgiN2KRpgvH/aXME8/dqc4Xr16avK7T4Ouc6MmfXWrvGWvuNtXYKMB5IxN2V4UQcxD2fp7YjIwt59ZwDoKiskqkLd3DeE9nMWneiA8oiIi1Wk/OuiIhTHF+gZq3Nxz2VodcJvnQV0N0YE1/reH/cO61sPNaLK12Wkooqbn91Kdtyi0/w0SIiLVKz8q6IiBMcL3aNMR2AfsCmE3zpB0A07ubUR+4VhbtJ9efW2rLG3KSiysXz2VtO8NEiIi2ST/KuiEh9yiqreHPxDp/fN6Bzdo0x04ClwAqgEOgD/BR327HHalx3ueePRxaxnW+MOQAcsNbOAbDWLjPGvAE8YYyJxr394I9x74J1bWNjqnRZpi3bxQMTBzTraxMRCXXHy72+yrsiIvW5793veHfpruNfeIICvUBtAXAl8HPcbWp24N6R6uFai9PeqvW6Zzwf5wBn1Tj+A+Ah4EEgBff+2ued6FauxeWVJ3K5iEi4akzu9UneFRGp7d7z+zFvQw7bfHxfx3dQC6QjPR+7/vpDr+OJsVGs/OO5jsR0RLCuqAzWuECxNUWwxgUhE5vaMpyAIzk3GP+fNkcwf682V7h+bfq6goe1ltnrDvDgR6v5y2WDyOrm3bHw6405nNG73ZFrQ7rPbtCIijBMGpLpdBhB+40arHGBYmuKYI0Lgjs2kZrC+Xs1XL82fV3OK6908dmqvbw8fyuLtrrbdf/feyv58K5RREUeXUJ2eq+2Pn92iy92oyMjuGV0d6fDEBEREQk75ZUuPli+mydnbmBbrvfmi2v3HuLLNfs4b0C6X2No0cVubFQEz1w3NGg2lhAREREJB+49DbbzfPYW9haWep0zBvp1TOIPF/dnRI80v8fSoovdUb3aMrZve6fDEBEREQkbL329lce/WE9BSUWdc5cOzeTHZ/akd4fWAYunRRe72RtyOFhcTpuEGKdDEREREQkLEQavQrdd61huPL0b147oQkp84GsuxzeVcMLgTskAlFe5mLbM9/3cRERERFqCPQUldY5dPbwLfTok0i0tnj9PGkj2r8Zyx9hejhS60EJbj722YBv3T/sOgD4dEvnsnjHV7TtERGpT67GmCdfWYyICW3KKeeSTtXy6ai8v3zScMX3aeZ3fnV9C+9axXp0WGsvXObdFjuxOOCWD+JhIANbvK2LZjnyHIxIREREJfpsPFPGzN79l7KOz+XTVXgBu+t8iPvf8+YiMlFZNKnT9oUXO2U2MjeKiQem8uXgnAG8s3MHQLm0cjkpEREQkOG3cf4inZm7kg+W7cdV6syarWyonZyY7E1gjtMhpDNZalm4/yKXPfA1AfEwkC39zNomxLbL2F5Hj0DSGptE0BpHQt/lAEY9/sZ6PvttD7R/lUzqn8McJJzO4c4pPn+nrnNtiq7shnVPo0yGR9fuKOFxexQfLd3P18C5OhyUiIiISFB7/fB3PzN5EZa2h3DF92nHHWT0Z3j01JNY8BcdkCgcYY5icdbS4nbpwu4PRiIiIiASX1IQYr0J3XL/2TLv9dF6+aTgjeqSFRKELLbjYBZg0JJMYz+Tp5TsLWL270OGIRERERAJv84Einpu72Wva0TUjutItLZ6RPVJ5/44z+O+NWQwJwTVOLbrYbZMQw3kDOlZ//sYije6KiIhIy7Ej7zAPfriacY/NYeba/V6jtTFREbz949OZcutIn8/LDaQWXewCTM7qXP3nact2UVpR5WA0IiIiIv5XUFLBQx+tZuyjs3l+3hbAPRe3qtb83LaJsSEzXaEhLb7YHdkjjS6p8QAUllbyyco9DkckIiIi4h+FpRU8PWsjY/46i+eyt3jNyV27t5AtOUUORucfLbYbwxEREYarsjrzt8/WATB14Q4mDenkcFQiIiIivlNaUcWTMzbwyoJtHCqt9Do3uFMyVwzrzOSszkGzEYQvtfhiF+CKUzvx+BfrqXJZvtmSx+YDRfRol+h0WCIiIiLNlr3hANe/sLDO8cyUVtx2Zg+uG9GViIjQnqpwLOFXvjdB+6Q4xvVrX/35G4t3OBiNiIiIiO9kpLTy+rx72wT+dvkgZv/yLL5/WrewLnRBxW61mgvV3lmyk/JKl4PRiIiIiJy4GWv2sXavdyvVnu0SyfQUvPee348ZPzuTK4Z1JjoMpyzUR9MYPM7s046OSXHsLSwlp6icmWv3cd6AdKfDEhERETmmiioXM9bs4+9fbGDdvkP8cEwP7r8gyeual27KIjMlnlYxkQ5F6ZyWUdI3QlRkBFcMO7owbeoiTWUQERGR4GStZdXuAv4wfRUj/jyDH726lHX7DgHu0d3aerVv3SILXdDIrpcrh3Xmn7M2Yi3MWX+AXfkl1cP+IiIiIk6rclneWbKTX72zosFrhnRpQ2WVKyw7KzSFit0aOqfGM6pXW7I35GAtvLV4B/ec3cfpsERERET4cvU+Hv18HWv3HqpzrmNSHBNOyeD6kV3p7Nk/QNxU7NZyVVZnsjfkAPDmoh3cNa43kWG+SlFERESC31ebcuoUulnd2nD3+D6c3jMt7LsqNJWK3Vq+178DqQkx5BWXs7uglOwNBzirb/vjv1BERETERw6VVhAdGUFc9NF5tteO6MqLX20lPiaSSUMyuezUTgzt0sbBKEODJnPUEhsVyaVDMqs/n7pQC9VEREQkMIrLKnn883UMf2gGn67c63WuV/tEHrlsIHN/NZaHJg1UodtIKnbrMXn40Z67X67Zx4FDZQ5GIyIiIuGusLSCv3+xnlGPzOTJmRspqahi2rJdda67KqsLbRNjHYgwdKnYrUev9q0Z1tX921Kly/LO0p0ORyQiIiLhqLSiiuezNzP6kVn8Y8YGDh6uqD63r7BUm1z5gObsNmDy8C4s3nYQgDcW7eC2MT0wRhO/RUREpPkqq1z896stPDt3CzlF3u8gt28dyy/P7ctlQztp0ZkPGGut0zEEjDHGgrsR8/EcLq9kxEMzOFRWCcDUH45kZI80/wYoIkHpyC+61lr9q3MCTiTnirQkX6zex6OfraveBOKI9q1juWt8by4bmkl8TMsdj/R1ztU0hgbEx0RxyZCM6s/f0I5qIiIi4gO5RWVehW7bxBjuO78fM39xFteP7NqiC11/ULF7DJOzulT/+ePv9lBQYx6NiIiIyPFs3F93A4jxJ3XAGIiLjuC2MT2Y88ux3HZmTxJjVeT6g/5Wj2FAZjInZySxanchZZUu3vt2Fzec3s3psERERCTILdt+kCdnbGDWugN1pkK2ax3Lv649laxubUhTZwW/08jucUwefnR0d8rC7Zp7JiIiIg36elMOlz7zFZOe+ZpZ6w4A8NTMDXWuO29ARxW6AaJi9zguOSWDuGj3X9PavYdYsbPA4YhEREQk2Ow8eJjL/vU11zz3DUu351cfNwZSE2Ipq6xyMLqWTdMYjiMpLpoLB2ZU99qdumgHgzunOByViIiIBION+w/xv6+38tbinZTV6IkbHWm4cGA6d47rRa/2rR2MUNR6rBEWbc3jin/PByAhJpKFvzmbBE0iF2kx1HqsadR6TMLdP77cwBMz1lP7W/ySUzK47/yT6Jgc50xgIU6txxwwrGsberZLAKC4vIqPVuxxOCIRERFxmsV6Fbr905N4+0en8Y/JQ1ToBhEVu41gjPFqQzZl0XYHoxEREZFA25ZbXOfY2Sd1AGBMn3a8evMIPvrJKIZ1Sw10aHIcmsbQSLlFZYx8eAYVVe7XfnbPGPp21BwckZZA0xiaRtMYJNQVl1XywfLdTFm4neU7C/jip2Po3eHov/3WWpZsO6gC18c0jcEhaYmxnNO/Y/XnUzW6KyIiEpYOHCrjgQ9XM/yhL7n33e9Y7unE9OzczV7XGWNU6IYAFbsn4KqsztV/nrZsF6UVaiMiIiISLnKLynj44zWM/utMXpi3heLyo//Ox0RGYNE7FaFILQVOwKhebclMacWu/BLyD1fw+ep9TBic4XRYIiIi0gz5h8t5LnszL361lcPl3gNZPdslcM2Irlw6JJM2CTEORSjNoWL3BEREGK7K6szjX6wHYOrC7Sp2RUREQtiOvMNc8I9sDpVVeh3vn57Ez77Xh/Enta+eQyqhSdMYTtAVwzoR4fme/3pTbr2rM0VERCQ0dGrTin7pRxed9e3Qmn9fN5QP7xrF2f07qNANAyp2T1B6civO6tu++vM3Fu1wMBoRERFpjCqXZdba/Tz44Wqv48YYfn5OX3q2S+Cpq4fwyd2jOW9AOhERKnLDhVqPNcHnq/byw1eWANCudSzz7x1HVKR+bxAJV2o91jRqPSbBoLiskv/M2cQ7S3exK7+ECAPLfncOya2iva5zuawK3CCh1mNBYGy/9rRrHQu425PMXLvf4YhERESkpkOlFTw7dxNjH53NkzM3siu/BACXhfmbcupcr0I3fGmBWhNER0ZwxamdeGb2JsA9leGckzse51UiIiLibweLy5m6aAf/nrOJgpIKr3Mp8dFcOqST18YQEv5U7DbRlcM6Vxe7s9btZ09BCenJrRyOSkREpGUqLK3g5a+38q/Zm7z64wK0TYzlnrN7c+WwzsRE6U3tlkbFbhN1a5vAaT3SmL85F5eFtxfv5K7xvZ0OS0REpEWau/4Aj36+3utYl9R4bj+rJxOHZBIXHelQZOI0/XrTDJOHH91R7Y3FO3C5tAhDRETECRcMSKd/ehIAGclx/OXSgXz5szOZPLyLCt0WTiO7zXDuyR1JiY8m/3AFOw+W8NWmHEb3bud0WCIiImFr8dY8/jlrI8O6tuHOcUffUY2IMPzu4v5szzvMxFMyNV1BqgX8O8EYc64xZqYxZq8xpswYs9MY86Yxpn+t6zobY942xhQYYwqNMe8aY7rUuqabMcY28F+Kv7+WuOhIJg3JrP58qnruikiYM8ac1UDOzXc6Nglf1loWb83jqv/M5/J/z2f2ugP8Z85m8g+Xe103skea5uVKHU6M7KYCS4BngANAF+BeYIExZqC1dpsxJh6YCZQBNwAWeBCYZYwZZK2tvW3Zw8D0WscO+fFrqDY5qwsvfrUVcPffzS0qIy0xNhCPFhFx0k+ARTU+r2zoQpGmKq908fo325iycAfr9nn/s36orJL3v93NDad3cyY4CRkBL3attVOAKTWPGWMWAmuBy4HHgFuBHkBfa+1GzzUrgA3AbcDjtW672Vq7wM+h16tvx9YM6ZLCsu35VFRZpi3bxS2jezgRiohIIK1xKu9K+KuocvH2kp38/Yv17D9U5nUuKsIwaUgmPzijO/0zkhyKUEJJsMzZzfV8PDIyMAFYcKTQBbDWbjHGfAVcQt1i11GTszqzbLv7HbwpC7dz86ju2ktbRESkCdbuLeSO15ay6YD3m7gxURFcMjiD28f2onvbBIeik1Dk2KQWY0ykMSbGGNMb+A+wl6MjvicDK+t52Sqgfz3HHzbGVHrm9043xgw8zrMb/K8pLhqUQUKMe6XnpgPFLNl2sEn3ERHn+DovtACvGWOqjDG5xpjXa6+pqEl/t3Ii4qIivQrdNvHRXDa0E/N+PZa/XTFYhW6YCGRecHJk9xvgVM+fNwLjrLVH9t1NBeqrGPOANjU+L8NdKH+Oe/5vP+B+4GtjzHBr7Rp/BF5bQmwUE07JYMpC9wK1KQt3MKxbaiAeLSISaAW4p5vNAQqBIbjz7nxjzJAaeVzkuAoOV7Az/zB9O7QmKtI9/tY1LZ5ObVpxsLicO8b14sbTuxEfEyxvREsoMtY60xvWGHMSkIR7bu4vgA7AKGvtVmNMOfC4tfbeWq95ELjXWtvgd70xpjPuEeDp1trrap2z4F7V6Wvf7shn4tNfARAXHcHC35xNUly0z58jIoF3ZKTBWquhyHoYY4YCC4G/WGt/W+O433KuhLb8w+X8Y8YGXvxqK49eMZjLT+3kdf7rTTn0ap9I+9ZxDkUoTvJ1znVsGoO1do219hvPgrXxQCLurgzgHtVtU8/LGhrxrXnfHcA8IMuH4R7X4E7J9Ovo3mu7tMLF+9/uDuTjRUQcY61dCqwnwHlXQk9pRRXPzd3MmL/Oqu5k9Njn6yiptb3v6T3bqtAVnwmKRnTW2nzcUxl6eQ6twj1vt7b+wOrG3tYHoTWaMYbJWTV2VFu0PZCPFxEJBhrClXoVlVXy7NxNjH10Ng99vIbC0qOd6vYUlLJoa56D0Um4C4pi1xjTAfd8202eQ9OBkcaYHjWu6QacQd1+urXv1QUYhfsttYCaNKRTdSPrlbsKWbmrINAhiIgEnDFmGNAXB/KuBLfSiiqez97M6Edm8ueP17KnoLT6XNe0eJ65dihbHr6AMX20+6j4T8Dn7BpjpgFLgRW4Fzf0AX4KdASGW2vXG2MSgOVACfBb3KMFDwCtgUHW2iLPvR7DXbDPx71ArS9wH5AMjLDWrqv1bL/PH7tn6jLe80xhuG5kFx6ceMzGECISAjRn9yhjzGvAFtx5PB/3ArX7gMPAUGttTo1rNWe3BXv/21389dN17Mov8TqemhDD3eN7c/XwLtrpTOrl65zrxPLGBcCVwM+BGGAHMBt42Fq7FcBaW2yMGQf8HXgFMMAM4J4jha7HKuDHwI245/zm4t557Y+1C91AmTy8S3Wx+/6y3fzmgv608rQlExEJAyuBq4G7gHjcbSPfBX5fs9AVKSmv8ip0O6e24vazejFhcAYJsequIIHjWDcGJwRilMFay9hHZ7M19zBAvatMRSS0aGS3aTSy23KUVVZRUWVJrFHEulyWy//9NVtyirl7fG+uGdFVI7nSKGHTjSFcGWO4Kutob3UtVBMRkXB1qLSC37+/kr6//ZT/zNnkdS4iwvD0tUPJ/vU4bjyjuwpdcYy+8/zgslMziYpw/zKyaOtBNu4/5HBEIiIivrN+3yHumrKMoQ98wUvztwHw0Yo9dUbx05NbeY32ijhBxa4ftG8dx/iT2ld//saiHQ5GIyIi0nzWWpZsO8jtry3hnL/P5YPlu6moOlrc5hSVcaCozMEIReqnObt+Mmvdfn7w4iLAvfJ0/n3jiI3SQjWRUKQ5u02jObvh44vV+3jk07Vs3F9U51zPdgncPKoHVw/vXP2zItIc4dCNoUUY07sdGclx7C4oJa+4nC9X7+fCQelOhyUiInJCfvHWct5esrPO8TN6pXH3+D4M757qQFQijadpDH4SGWG4YtjRHdWmaqGaiIiEoIsHZ1T/OSEmkgmDM5h2++m8dstIFboSEjSNwY92HjzM6L/O4sjjsn81ls6p8QF5toj4jqYxNI2mMYSWfYWlvLV4B13SEphQo8AF94ZJpRUu/nLZQFLiYxyKUFoKTWMIIZ3axDOmdzvmrD8AwJuLd/Dzc/o6HJWIiMhRmw8U8fSsTUxfvouKKsv3+nfwKnYrq1zcc3YfuqbFa06uhCQVu342OatzdbH71uKd3D2+N1GRmj0iIiLOKa908ebiHbz2zXbW7Cn0OvftjnystdWFbVRkBN3aJjgRpohPqNj1s/EndSAtIYbc4nL2FpYyd8MBxvXr4HRYIiLSAhWVVfLK/G28+NUW9h+q2yZscKdkLh/WGWtBg7gSLlTs+llMVASXn9qJ/8zdDMCUhTtU7IqISEBZa3lh3hYe/GhNnXPGuDsI/WR8b07t2saB6ET8S8VuAFyZ1bm62J25dj/7C0tpnxTncFQiItJSuCzVU+qOSIqL4tbRPbhmRBfSEmMdikzE/zR5NAB6tkusbs9S5bK8VU+/QhEREX+JjDA8dfUQung6At13fj8W3D+eu8b3VqErYU/FboBMzjrac/fNxTtwudSKR0REfKvKZZm+fDeX/+tr9haUep1LiY/hX9cNZfFvz+a2M3sSH6M3d6VlULEbIBcMTKd1nDuxbMs9zILNuQ5HJCIi4cLlssxYs4/xj83mJ1OWsXjbQX4/fWWd607OSKatRnKlhdGvdQESFx3JpCGZvDx/GwBTF+3g9F5tHY5KRERCWW5RGa8s2MarC7aTU+TdXWHO+gNsyy2ma5rahknLpmI3gCZndakudj9duZeDxeW0SdBONCIicmIKDlfwwrzNPJu9mdIKl9c5Y+DG07vxwzE9SE9u5VCEIsFDxW4A9c9IYlCnZFbsLKC8ysW0Zbu4aVR3p8MSEZEQUeWy/GH6Kt5cvIOySu8iNzUhhvMHdOTu8b3V8UekBhW7AXZVVmdW7CwAYOqi7fzgjG7aflFERBolMsKQU1TmVej2T0/iR2f15PwBHYnWDp0idRhrW05XAGOMBXdzbaccKq1g+EMzKKmoAuDd209naBc18RYJZkd+IbXW+vU3U2NMLHAaMBLIAFoBOcA6YK61drM/n+9rwZBzQ1lllYtPVu4lq1sqHZOPjtRu3H+Ic/4+l34dk7hpVHcuHZJJRIQGTSR8+DrnamQ3wFrHRXPRoPTqXrtvLNyhYlekhTPG9ALuAa4FkgEXUACUAKlAHGCNMUuAZ4CXrbWuBm4nIa6kvIqPvtvD07M2siWnmFduHu5V7PZq35rP7hlDr/aJemdQpBH0focDJg/vUv3nD1bspqis0sFoRMRJxpingdVAFvAnz8c4a22atbaTtTYeSAcuBb4FHgdWGWNGOBWz+EdllYt/ztzAGY/M5BdvLWdLTjEA+wrL6lzbu0NrFboijaSRXQcM7ZJC7/aJbNhfxOHyKj5YvpuraxTAItKiZADDrbXfNnSBtXYf8D7wvjHmLuA2YDDwTWBCFH8qKa/inaU7eXbuZrbnHfY6lxQXRUqraMWvRloAACAASURBVIciEwkPmrPrkBfmbeGBD1cDMLhTMu/fOcrhiESkIYGasxtuginnBqNDpRVMXbiDp2dvJP9whde51rFR/Oisnlw7ogsp8WpRKS2Lr3Ouil2H5BWXM/LPMyivck+7+/gno+mfkeRwVCJSH38Wu8aY7wNvW2sPH/fiEBNMOTfYfLUxh+te+IbafzWtY6MYf1J7/jRxAElxGtGVlknFbjMEW+K9a8oyPli+G4DoSENllSUhNoqJQzK4dXQP7XojEiT8XOy6gCLgXeAla+0sXz/DKcGWc4OJtZY7pyzjoxV7AOic2oqbzujOFcM6kxirGYbSsqnYbYZgS7xPz9rI3z5bV+d4VIQhOjKCZ64byti+7R2ITERq8nOxOx74PjAJSAB2Aq/g7riw3tfPC6Rgy7lO2XnwMMu253Px4Ayv4/mHy7nxxUVcNCid75/WjZgorRkXARW7zRJMiXdbbjHnPZFd3W+3Pq2iI/n0ntEa4RVxWCDm7Bpj4oHLgOuBcYABFgIvAVOttfn+era/BFPOdUJuURmPfLqWd5fuIiLCMO/XY2nfWjubiRyPr3Oufo10yHPZm6moOnabzIoqF89nbwlQRCLiJGvtYWvtK9bac4DOwH24R3qfAfYYY94yxlzkaJDSKLvzS/jzx2s482+zeXPxTipdlvJKFy/MUz4XcYJGdh0y4PefNaq/bmJsFCv/eG4AIhKRhjjZjcEYcwrwI+BWdwg2ZCZ0BlPODYSVuwp4PnszH67YQ6XL+2se2SOVn4zrzem92joUnUjo0A5qYaK4kRtJFJdrwwmRlsoYMw73tIZLcU9r2OhsRFKbtZY56w/wXPZmvtqYW+d8z3YJ3HN2nzrzdUUkcFTsOiQhNqrRO6e9u3QnEwZnEBWpWSci4c4YcxLuBWvXAplAITAFd6eG+U7GJnUt3JLHjS8uqnN8RPdUfjimB2P7ticiQu2ZRZykaQwO+e173zF14Y46b3U1pGtaPLef1ZNJQzppxa5IgPl7GoMxph1wDe5R3CGAC/gc9+K09621dfeLDQHBlHP9xVrL+f/IZu3eQ0QYuGBgOreO7sHgzilOhyYSstSNoRmCKfE2phtDfTJTWvGjs3py5bBOxEZF+ik6EanJz63HPgTOwf1O20rcBe5r1tq9vn5WoAVTzm2u4rJKXp6/jcGdkuvMu/1wxW6WbDvITWd0p3NqvEMRioQPFbvNEGyJd9a6/dz+6lIqqlxeI7xH+uw+evkgNhwo4r/ztlBY6j3loUNSLLeN6cnVw7vQKkZFr4g/+bnY3Q+8jnuawjJf399JwZZzm+JgcTlvLN7Bs3M3k1dcTq/2iXxy92iiNa1MxG9U7DZDMCbebbnFPJ+9hWnLdlFcXklCTBSThmRyy+ju1f11D5VW8PL8bbwwbwt5xeVer2+bGMOto3tw3ciuJGjXHRG/8HOxG2WtDcuVqMGYcxvDWsuq3YX8Z+7m6l0ua3r40oFcPbyLA5GJtAwqdpshVBPvEYfLK3ltwXb+M3czOUXeU/jaxEdz86jufP/0btpPXcTHAtl6zBiTCfwcGAOkAhOstSuNMfcA86213/g7Bl8JxZz7zeZcfvjKEgpKKuqcy0xpxV3jenHpUK2dEPEnFbvNEIqJtz6lFVVMXbidf8/ZzN7CUq9zSXFR3HhGd246oxsp8TEORSgSXgJV7BpjTgaygSpgPnAhkGWtXWqM+TvQwVp7jT9j8KVQyrkFhyt49PN1vLJgW51zKnJFAkvFbjOEUuJtjLLKKt5espN/zd7EzoMlXucSY6O4/rSu3DKqO2mJsQ5FKBIeAljsfgq0Bs4FSoFyYJin2L0CeMRa28OfMfhSKOXcXfklnPGXmV7Hzj25A+P6teeSUzKJi9baCJFAUbHbDKGUeE9ERZWL95bt4pnZm9iSU+x1rlV0JNeO6MIPx/SgfZL2ZBdpigAWu0XA1dbaD4wxkUAFR4vdMcCn1tqQWe4frDm3sspFcVkVyfHeU74ueforlu/IZ1CnZH53UX+GdUt1KEKRlk3FbjMEa+L1lcoqFx99t4enZm5k4/4ir3MxURFcndWZ287sSUZKK4ciFAlNASx2C4FrGyh2LwWes9am+TMGXwq2nFvlsny4Yjd3T/2WS07J4B+Th3id//i7PVRUuZgwOKP6/7mIBJ6K3WYItsTrLy6X5dNVe3lyxgbW7j3kdS460nD5qZ25/aye6gcp0kgBLHa/BAqttZfWU+xOBeKttRP8GYMvBUvOPVhczv++3spbi3ewu8C9ziElPppv7h+vfuUiQUjFbjMES+INFJfLMmPtfp6auYEVOwu8zkVGGCYNyeSOsb3o3jbBoQhFQkMAi90zgS+BWbh7774A3AecDEwGxqgbQ+Mt236QZ2Zv4ovV++o9//aPTtNUBZEg5Fixa4yJBU4DRgIZQCsgB1gHzLXWbvZFQP7kdOJ1irWWOesP8NTMjSzZdtDrXISBiwdncOfYXvTu0NqhCEWCW4Bbj10IPAH0rHF4K3CHtfYTfz/fl5zIuev2HuKtxTuYtW4/mw4U1znfOi6Kiwal88MxPfWLvkiQCnixa4zpBdwDXAsk496zvQAowd0DMg6wwBLgGeBla63LF8H5Wkstdo+w1jJ/Uy5PztzAgs15XueMgfMHdOTOsb3pn5HEttxinsvezHvLdlNcVklCbBQTh2Rw6+ge1ZtdiLQUgSx2azyzF9AeyLXWrgvUc33JiZxbVFbJhKfmsbnWYt1ObVpx59heTBqaqakLIkEuoMWuMeZp4FZgGTAVmAssr7nbjzGmA+7R3guBy4F9wI3B+FZbSy92a1q0NY8nZ2wge0NOnXOndE5hzZ5Cqly23m2Mn7luKGP7tg9kuCKOCuA0hoHW2u+Ocf4qa+0b/ozBl/ydc0srqjAGr+J15a4CLnpqHuDuRnP+gI784IzuDOyU7JcYRMT3Al3sTgP+aK39tpHBxQK3AaXW2md9EaAvqdita9n2g/xz5kZmrN3f6Ne0io7k03tGa4RXWowAFru7gNOstdvrOXcl8Kq1NmR2i/FXzq2ocvHPmRt5af5W/nDxyUwckll97s1FO/hk5R4uO7UT4/t1oFWMRnFFQo0WqDWDit2GrdxVwD9nbuTTVXuPe21UhOHq4V14YOKAAEQm4rwAFrsfAL2B0621eTWOX457wdrj1tp7/RmDL/k65+4rLGXqwh3896st1dv5ntmnHS/dNNwn9xeR4ODrnHvMfQ+NMac29kbGmH82PxxxyoDMZP59/anEN2IUpNJlmbZsVwCiEmlxrgTygI+NMfEAnv66rwNPhlKh60s5RWX8c+YGxj82h79/ub660AXYlltMeWVQLhMRkSBxvGkMOcBZ1tqVx7yJMS/gnqcb1O8XaWT3+Lrf+xGN/du5dkQXJgzOIKtbKhERasAu4SvA3RhSga9wd2D4L/Aq8Iy19qf+fravNSfnWmtZuCWP95fv5qMVe7wKXHBviX7H2F78cEwPIpV/RMJKoOfsLgPSgTPrWw1s3NG8jLtTw/9Zax/yRVD+omL3+Ab8/jOKyiqPf2EN6clxXDQonQmDMxmQmaSdhyTsBLobgzGmM/A17jaPT1trfxKI5/pac3LuHa8v5aMVe+oc75YWzx1je3Hx4AziooN6fEVEmijQxW4aMAdIwV3wbqpxLhL3W2tXAL+y1j7qi4D8ScXu8f32ve+YunCHVxeGE9G9bQIXD85gwuAMerVP9HF0Is7wZ7FrjPlTA6dOBkYD/4HqN1ystfb3vo7BX5qTc5/P3syDH62p/jwzpRV3juvF5ad2IjrymDPwRCTEOdFntz3ulmNxuHfv2W6MiQbeAiYAd1trn/JFMP6mYvf4tuUWc94T2ZRUVDV4TavoSP5y2UAWbsnj4+/2cPBwRb3X9U9PYsIpGVw8OIPMlFb+ClnE7/xc7J7IhFMb7NPFampMzrXWsmp3IQMyvVuDFZdVMuavsxjTpx3nntyBcf06EBOlIlekJXCkG4MxJhP3CC/A+cCTwDnAj4OxxVhDVOw2zqx1+7n91aVUVLmO22e3osrFvI05fPDtbj5btZfi8vqL5GFd2zDhlAwuGJhO28TYgHwdIr7ixKYS4eB4OXfZ9oP85ZO1fLMlj/fvOIPBnVO8zheVVZIYG+X/QEUkqDi5XXAX3AVvJmCAm621L5/wA405F/g10B9oAxzAPTftD9ba1TWu6wz8Hfie53lfAvfU7j9pjGkD/A2YiHsL4/nAT+trzK5it/G25RbzfPYWpi3bRXF5JQkxUUwakskto7s32F+3tKKKGWv2M335LmatO1DvCunICMPpPdOYMDiDcwd0JCku2t9fikizqdj1dgL5uU7OrXJZvlyzj5fnb+WrjbnVx0/rkcbrt47QnH8RCfic3ZtqHeoJ3Ad8ArxT+3pr7X+P+0BjrgaGAt/gLnS7APcCnYGB1tptnpY7y4Ey4Le456s9CMQDg6y1xZ57GSAb6Ab8Ejjoie9k4BRr7c5az1axGyCFpRV8vmof05fv5quNOVTVMwc4JiqCsX3bMWFwJuNPaq/FJhK0/DyNIc5aWxqo1zVXY/Oz59rqnFtZ5eLpWZv4+5fr69zzSO/u/7uov6YqiEjAi92AzCUzxvQF1gK/sNY+Zoy5G3gc6Gut3ei5pjuwAfdiuMc9xy4B3gPGWWtneY4lA1tw7zT0k1rPUbHrgJyiMj75bg/Tl+9m0daD9V6TEBPJOSd3ZMLgDEb1bqsFKBJU/Fzs7gUeAV601uY34vrTcQ8QLLLWPuDreBrx/EblZ89xC/D0rA28tmA7u/JLvO4VYeDiwRncPb43PdppQauIuAW62O16Ijez1m5rUhDGtMU9ynuPtfYfxpgZQJy19oxa183xPOdMz+cvAOdZazNrXfcS7v7AXWsdV7HrsF35JXy4fDfTl+9m1e7Ceq9pEx/N+QPTmTA4g+H19PDdllvMc9mbeW/ZborLKkmIjWLikAxuHd1DWxiLX/i52L0U+DPQFfgU97tVy3HnxDLc0716AMOBi3C/G/Yi8Dtr7T5fx9OIeBuVnz3HLEDXX39Y5z5XnNqJH5/VU0WuiNQRNtsFe1qXReJO8H8BTgcGW2v3e0Y63rfW3lbrNc8AV1hr23k+XwAUWGvPrXXdr3CPlLS21hbVOH7cL1aFcOBs3F/EB8t388Hy3WzOKa73mo5Jnh6+p2QwMDOZ2esPNHrxnMiJaMxcUX/N2fXkw4nAzcBZuLvf1ExGBtgGvAE8a63d7I84GqOx+dlzzKvYjSo5RP7KGRyc/yauEvcvu8q5Ii1TIHOuk8tcvwGObEe8EfdUhP2ez1Nxz7+tLQ/3KAc1rtvawHV4ri2q57wEgV7tE/np9/pwz9m9WbW7kOmewndPwdFpiHsLS3l+3haen7eFzJQ49hWW1dsDuNJlqXRVcfurS/n0ntEa4ZWQYq2twr0O4h1jTAxwCu4NJeKAXGCttXaHgyHW1Nj8XC2tJJ/r137OmTuXYSrKmJIYx+MlhZzY9jUiIk1zzGLXGDMd+L21dlljbmaMiQNuBw5ba/99nMuvB5Jwvz33C+ALY8woa+3WxjyrOTSSEFyMMQzITGZAZjL3ntePxdsOMn35Lj7+bi95xeXV1+3KP/5anIoqF89nb+GBiQP8GbKEoeNM6QpkHOXAwoA9MABemf8kptAzbSkightTU7nl5JNJ+8EPcJWUENFKfbhFWppA5tzjrQLaCiwwxnxjjPmJMWaoMcarQDbGZBhjJnrmz+7B/Tbc0uM92Fq7xlr7jbV2CjAeSMS96ALcowb1jRDUHlE41nVQ/+iDBLGICMPw7qk8OHEg39w/nv/9IItLh2Y2utdmpcsybdkuP0cp0qI1Nj9X6ztzBu1/8XOvY5V79rDvz39mw5gzOfD007hKSup7qYhIszVmB7WewD3AtUAy7nlkhbgXTqQAMbjnky0E/oW7C0LD2281/JzFQL619mxjzEwgxlo7qtY1sz0xH1mg9l/gHGttp1rX/Q8YqwVq4aO0oop+//dpo6+/e3xvRvduy+DOKersIM2mPrtHNTY/e4555VxXSQn7HnmEwk8+xVVQ4HXfiKQk0m66idTrryMiQdOQRFoyJzeViAFOA0ZQay4ZMLepnRg89+4AbAJes9beZoy5B3gU6HNkIYYxphvu1jb3Wmsf8xybCEzD3XlhjudYEu7WY69ba++q9RwVuyFswO8/o6jsxGb5JcZGMbJHGqN6pTGqdzt6tktQ03o5YSp2j2psfvYcrzfnVublUTBtGrkv/o+qnByvc5Ft2pB28020ufpqFb0iLVTId2MwxkzDPc1hBe4R4j7AT4GOwHBr7XpjTALu1jslHG1a/gDQGnfT8iLPvSKAebg3pKi5qcQg3J0dvBZ0qNgNbb997zumLtxR7wK1xspIjuOMXm0Z1bsto3q1JU1bF0sjqNg9qrH52XPtMXOuq6yMvP/+l7zXXq9T9JpWrUgYOZLUG24gYeQIv3wtIhKcwqHY/TVwJe7d2GKAHcBs4OGai9M82xPX3I5yBu4+vFtr3S8V9yjDRNyjzfOBn1lrl9fzbBW7IWxbbjHnPZFNSUXDs2TioiP41bn9WLf3EPM25tRpYl9b//QkRvd2F79Z3VK1i5vUS8WutxPIz43Kua6yMgrefZfc556nYvfuOueTLjif9r/6FdEdO/rmCxCRoBbyxa6TVOyGvlnr9je6z661li05xczbmEP2hhwWbMrl0DGmQcRGRZDVLbV61Ld/elKdDS2kZQpUseuZD3u7tXZtPef6AP+21o7zZwy+dKI515aXkz/tPXJfeIGK7dvrnE844wzSbr2V+BHDNR1JJIyp2G0GFbvhYVtuMc9nb2Hasl0Ul1eSEBPFpCGZ3DK6+zH761ZWuVi+M5/sDTnM25DDsh35VB1jSkRaQgyn92rLaM+0h4wUtUdqqQJY7LqAkdbaOq3HjDGnAgubui27E5qac21VFYe+nEHuf1+gdPmKOudjevQg5dJJpFx+OZEpKb4JVkSChordZlCxKzUdKq1gweY85m04QPbGHDYfqH8XtyN6tEvwFL7tGNkjldZx0XWu0VbG4SnAxe4Ia+2ies5dCTxnrU32Zwy+1Nyca63l0CefkPPc85StWVPfA2h7xx2k3ngjkYn6+RIJFyp2m0HFrhzLrvwSvtqQQ/bGHL7amOO1oUVtURGGUzqnMKp3W3eLs04pZG/M0VbGYcqfxa4x5gfADzyfnoF78e6hWpe1AgYAM6y1F/k6Bn/xZc4t27yZvFdeoeCdd7Hl3j+bEYmJpFx+OSlXXkFsjx7NfpaIOCugxa4xZjMwqb7FXqFIxa40lstlWb2nkHkb3VMeFm7No7zS1eD1CTGRlFRUcaxGEa2iI7WVcYjyc7F7A3Cj59MzgWW4O9XUVAasBh6x1u7zdQz+4o+cW5mTQ/7bb3PwzTep3L2nzvlWw04lYfhw2v7oR5iYGJ89V0QCJ9DFboPzx0KRil1pqtKKKhZtzWPeBvdit9V7atcixxcVYbh6eBdtZRyCAjiNYRbw4/oWqIUif+ZcW1FB/jvvkvfSS5Rv2VLnfGTbtqRedy0pV11FVJv6NnwTkWClYrcZVOyKr+QUlfGVZ9R33sYc9hSUNup1sVERvHv76fTrmESkOj2EDLUea5pA5FzrclGcnc3B16dQNGdO3RhiY0k6/3ySL51EwvDhfotDRHzHiWK33sUSoUjFrviDtZYe933MiXxXtY6LYljXNmR1T2V4t1QGdkomNipkFtm3OIEudo0xg4G+uHuHe7HWvhyIGHwh0Dn38NJl7H/8McrWrsNVVFTnfKvBg0m+9FKSL7pQu7OJBDEnit1PgZwGLzrKWmtv8EVQ/qJiV/ylKVsZ1xQbFcHgzimM6J5KVrdUhnZtQ2JslA8jlOYI4DSGFOAjYOSRQ56P1UmrJbQeay5bXk7hJ5+Q+9JLlK2u28XBxMW5pzhMnkxMp04BjU1Ejs+JYncv7sURx2OttUG9DFbFrvhLY7YyjjDQuU08hyuqOHDo2D9SkRGG/ulJZHVLZXj3VLK6tdHWxg4KYLH7DDAOuBnIBiYBBcBNwGnAZGvtEn/G4EtO51xrLYe/Wcj+xx+ndPVqqKz1C2l0NCkTLyHtlluI6drVkRhFpC7N2W0GpxOvhK/GbGV8pBtDl9R4tuUeZuHWPBZuyWPR1jy25R4+7jN6tkvwFL7uArhTm3hffglyDAEsdjcBfwReAyqArCPFrTHmX0CCtfb7/ozBl4Ip51bm5VEwbRo5//p33SkOkZEknnkm7e65m7g+fZwJUESqqdhthmBKvBJ+TmQr49r2FZayaGsei7bksXDrQdbuLeR436YZyXFk1Sh+e7VLbHB7Y2120TwBLHYPA+dYa+d5/nyBtXa259z3gKnW2jR/xuBLwZhzXWVlHHx9CvlvvUX55s3eJ40hPiuL9Af+pJFeEQep2G2GYEy8El6aupVxbQUlFSzZlsfCLQdZtDWPFTvzqag69vdtm/hohnVzL3jL6p7KyRlJREdGNKsIF7cAFrubgZ9Yaz80xqwBXrbWPuw592PgQRW7vmGrqjj8zTfkPPschxcsqHM+/rSRpD/wIDGdMh2ITqRlU7HbDMGceEWOpbSiim935FdPe1iy7SCHyxueMgEQHxNJv/TWrNhRcMy5xNrs4vgCWOy+DGy31v7WGHM/8DvgJaASuAGYbq29xp8x+FKo5NySb79l/xP/qLfojWrfnvSHHiRh1Kjq7wMR8S9tF9wMoZJ4RY6nssrF6j2FLNzinve7eNvBY25vfCza7OL4Aljs9gQyrLXZxpho4C/AVUA87s44d1lrc/0Zgy+FUs611lI0azb7H3uM8k2b6pxvNXQoaTffROLYsZiICAciFGk5VOw2QyglXpETYa1l04Gi6mkPC7fksSu/pNGvT4yNYuUfz/VjhKFNm0o0TSjmXGsthdOns/+xx6ncv7/O+diTTiLtlptJuuACjfSK+ImK3WYIxcQr0lS78ks44y8zG339LaO6c8GgdIZ0TtE/4rU4UewaYxKBNGC3tbYiUM/1pVDOua6yMgqmTyd/6huUrlpV53xURjppN91MypVXEBET40CEIuFLxW4zhHLiFWmKpmx2kZnSigsGduTCQRkM7pSswpfAFrvGmIuAPwGDPYeyrLVLjTHPAzOtta/7OwZfCZecW7FnD3mvvsrBV1/Dlnn3yI7q2JG0W24heeJEIhM1713EF1TsNkO4JF6RxmrMZhfHkpnSigsHpXPhwHQGteDCN4BzdicC7wAzgM+BvwLDPMXub4Ax1tqQmW8Sbjm3Yt9+cl94noJp7+E6dMjrXERCAsmXXELqjTcQ06WLQxGKhAcVu80QbolX5Hgau9nFHy85mcVb8/hs1T4KSup/x7xTm6OF78DMllX4BrDYXQYssdbeYoyJAso5WuxeAjxjrQ2ZXljhmnNdpaXkv/EGOc8+R1Wu93pBEx1N63POIe22H2qDCpEmUrHbDOGaeEWO5UT67JZXuvh6Uw4frdjDZ6v2Ulha/xSIzqmtuHBgBhcNSufkjKSwL3wDWOyWAhdba78wxkTi3kXtSLE7BvjcWhvnzxh8KdxzrqukhPy33+HglCl1N6gAki44nzZXX02rYcPC/mdExJdU7DZDuCdekYY0ZbOL8koXX23M4aPv3IXvoQYK365p8Vww0D3iG66FbwCL3f3A3dbaKfUUu98HHrLWdvZnDL7UUnKutZbDCxZw4B9PUvLtt3XOx3TtStvbf0zSRRdhIiMdiFAktKjYbYaWknhFfK2ssoqvNubw4Yo9fLFqH4caWPTWLS2eCwelc8HAdPqnh0/hG8Bi9zVgIDAGOIS72D0VWA1kA99aa3/ozxh8qaXlXGsthxct4sCTT1KyeEmd87H9+tHxN/cTn5XlQHQioUPFbjO0tMQr4g9llVXM2+Ce6vDF6oYL3+5tE7hwYDoXDkqnX8fWIV34BrDY7QYsBCzwMfB94G1gEJCMe5R3tz9j8KWWnHNL164l7+VXKJg+HSq9f0YSzzyTNtdfT8IZp4f0z4WIv6jYbYaWnHhF/KG0oorsDTl8/J278G2ozVmPdkcL374dvAvfbbnFPJe9mfeW7aa4rJKE2CgmDsng1tE9gmYL4wC3HusE/BE4F2gP5OLePe131tod/n6+LynnHu3gkP/Gm3XalsUNHkTH//sdrQac7FB0IsFJxW4zKPGK+E9pRRVz1x/go+/28OXqfRSX198Bome7BC4c5F7ctiu/pNGL55ykHdSaRjn3qIrdu9n/xBMUTv+gzrlWgweTfOmlJF1wPpGtWzsQnUhwUbHbDEq8IoFRWlHFnPUH+GjFHr5cs4/DDRS+Bvf79Q1pFR3Jp/eMdnyEV8Vu0yjn1lW2aRMHX3uN/LfexlZ4t/kzMTEkXzKBtnfeSXSHDg5FKOI8FbvNoMQrEnilFVXMXuce8Z1xjMK3PlERhquHd+GBiQP8GOHx+bPYNcb87gQut9baB3wdg78o5zasfPt2dt9/f70L2QBSb7yRdnfdSURCcEzlEQkkFbvNoMQr4qyS8ipmr9vPR9/t4cMVexr1msTYKFb+0dlNw/xc7LrqOWxxD3zXOW6tDZneVcq5x1eZl0fhhx+RP20aZWvWeJ2L6tCB1O9fT9LFFxPd3vnpPCKBomK3GZR4RYJH93s/OuYUhiOMgS0PX+j3eI4dg1+L3drFaxRQAowAlta+3lrb+KFxhynnNp61lkOffMKun/283vNtrrmGdvfcTWRSUoAjEwk8X+fcCF/cRETkRCXERjXuupjGXReqrLVVNf8DjrS0qKp9LpQKXTkxxhiSLriAviuW0+H++4ls19br/MHXX2f98BFsv+lmynfudChKkdCkYldEHDFxSAZREcf+pT0qwjBpSGaAIhJxXkRMDKnfv55en39O+p//TEzXrl7ni7/+mk1nf4/d/wdi/gAAIABJREFUv76Xytxch6IUCS2axiAijtiWW8x5T2RTUtHwYGVL7MZQe5tgfz/Pn5Rzm89WVpL3yqvkPPUUrsOHvc6Z6GgSxowh+cILSDz7bCJiYhyKUsS3NGe3GZR4RYLLrHX76+2zCxAZYXj+hmEtrs+uil2pjy0vJ+fZ5yh4/30qdtTdWyQiOZm0m2+mzTXXEJmoDg4S2lTsNoMSr0jw2ZZbzPPZW5i2bBfFZZXVi9ZaRUew8Ddn0zou2tH4wO8L1HrUOhQJrAMuAVbVvt5au9nXMfiLcq5/HJo5k5z//IfS5SvqPZ9w5hgyHnqIqLZt6z0vEuxU7DaDEq9IcLPWcu4Tc1m/rwiA31/cnx+c0d3hqALSeqx2Umpwvw21HpMjyjZtouD96eS/8w5VtebvRrRuTbuf3kPK5ZdreoOEHBW7zaDEKxL8Xl2wjd++txKA7m0TmPGzM4k4zkI2f/NzsXvDiVxvrX3J1zH4i3JuYNjycvJee528F1+kcv9+r3MRSUm0ueoq2lx7DdEdOzoUociJUbHbDEq8IsGvuKySkQ/P4FCpuwPX/36QxVkOz9vVdsFNo5wbWK7SUvL+9xK5zz5bZzEbQMrVk+l4//2YaOenBokci/rsikhYS4iN4sphnas/f+nrrc4FIxJCIuLiaPuj2+g1dy7tfvpTTGys1/n8KVPZcNZY9j/2GGUbNjgUpUjgaWRXRILOttxiznp0Nta6d1Cb9fOz6NbWuRXmGtltGuVcZ9nKSopmz2bfX/9Gxfbtdc7HZ2WR/sCfiOnWLfDBiRyDRnZFJOx1TUuobjlmLbw8f5vDEYmEHhMVReuzz6bnhx+QcvVkIuLjvc4fXrSITeedz66f/ZyiuXNxlZU5FKmIf2lkV0SC0pz1B7jhvwsBaB0bxYL7xzd6i2Ff08hu0yjnBhdXaSnFX88n99lnKfn22zrnI1NTaTN5Mqnfv57IlBQHIhRx0wK1ZlDiFQkdLpfl7MfnsDmnGIAHJg7g+pFdj/Mq/1Cx2zTKucHJWsuhL74g76WXKVmypM55Ex9PfNYw2t15F60GDnAgQmnpVOw2gxKvSGj531db+MMHqwHo3T6Rz386pjoJBpKK3f9v787Do6yuB45/TzZI2BJ2QQURq2JVUKu4AiKoYAUUtwpo3bFurYjy0ypaq1gV94W6UgEFq6K4obKIyqIoKlBBUVYBIZAFQiDb+f1x38BkMgnJzCTvZHI+z5Mn5L3vct4JOTm5c997w2M5N/bt+Pprst94k20zZlCSk1OuvcEhh9D2rjtJ69bNh+hMfWXFbgQs8RpTt2zbWUj3+2aQV1AMwMQrjuPEzrW/KpQVu+GxnFt3aFERudOnk/nkUxSsXFmuPe2442g+bCiNe/VCEuxxH1Oz7AE1Y0y90aRhMoOP3nf31y/bNGTG1AhJSqJZ//50eudt2t59d7n2HQsWsO4v1/HLH88m98Pp9geMqVOsZ9cYE9N+3ryd3g9/CkCCwKe39GK/5ml7OSq6rGc3PJZz67YdixaR9coEcqdPh+LiMm0pBxxAuwfGkHrEET5FZ+KZDWOIgCVeY+qmoS8s4LOfMgG4+pROjOp3aK1e34rd8FjOjQ+FGzaQNXEiWa++RkleXpm2hr//PRlDLqbZgAG+jKc38cmK3QhY4jWmbprxw29cPn4hAM1Sk5k/qjepKYm1dn0rdsNjOTe+FGVmsvHuu9n28Sfl2hLS0mjcuzdtbruVpBYtfIjOxBMbs2uMqXd6Htya/b2hCzn5hbz97a8+R2RM/ZPUsiX7PvEEnd6dRtrx3SHgQbWSHTvInTaNFaf2ZtOjj1KUleVjpMaUZT27xpg64fnPfuHe934A4JC2TfjgxpNr7W1T69kNj+Xc+FaUlcXWF19iy3PPlWtLSEujSZ/TyBg2jNTDDvMhOlOX2TCGCFjiNabuysl305DlF7oHZSZf1Z3jOtXO26VW7IbHcm79oAUFbBk/npy3plLwyy9lG0Vo2q8frW64npQO/iwKY+oeG8ZgjKmXmqUmM+io9ru/Hj9vlW+xGGP2kJQUWl55JZ2mvUP7sQ+T0vnAPY2q5L73Hr+c9Uc2jX3EhjcYX1jPrjGmzli+cRunPzoHgMQE4bORvWiXnlrj17We3fBYzq2ftKSEbR99xG9jHqBo48Zy7S2vvZaMIReT1Ly5D9GZuqBO9+yKyGAReUNEVotIvogsF5H7RaRJ0H5dReRDEdkuIrki8o6IdA5xPq3go2vt3ZUxprYc3LYJx3tDF4pLlIkLVvscUf0jIqsqyLsD/Y7NxAZJSKDpGWdw0OxZdJwymdQjjyzTnvn006zo2Yv1t40if8lSn6I09Umt9uyKyHxgDfA2sA7oBowGlgEnqGqJiBwEfAMsAcYAScBdQEugq6puCjifAi8D44Iu9b2q7ghxfetlMKaO+3DJRq6Z8DUAzRulMPe2U2mYXLPTkFnP7h4isgqXs0cHNS1X1aygfS3nGrSkhNxp09wwht9+K9fe6ITjaT1yJA0POcSH6EwsqtMPqIlIK1XdHLRtGDAe6K2qM0XkeWAw0FFVs7199gVWAI+r6siAYxX4p6reUcXrW+I1po4rKi6hx4Oz+TU7H4CHzjuyzJLCNcGK3T28YvdzVR1ShX0t55rdtKiInHffJWviJHYuXly2MTmZFpcMo8VVV5HYtKk/AZqYUaeHMQQXup6vvM+lT550B+aVFrrecetwPb2DajZCY0ysS0pMYEj3PU91j5+7yoopY+oASUoifeBADnh9Ch0nv0bT/v33NBYWsuX5F1jRpy+Zz46jKDPTv0BN3ImF2Rh6eJ9/8D4XAwUh9tsFHCgiDYO2DxeRXSKyQ0RmisjJNRWoMSY2XPiH/WiQ5NLX4l9z+GZN9l6OMFH2Ry/n7hKR+TZe11RX6pFH0v7hhzjg7amkduu2e3tJTg6bH32Un884k81PP22zN5io8LXYFZH2wD3AJ6q60Nu8HDhaRJID9msCHAYIkBFwignAtcBpwFVAC2CmiPTcy3Ur/DDGxL6MRikM6Npu99fj566K+JyWF6psGnA9cDpwMbATeEtEKhzWYK+tqUjDgw+mw6SJtHvwXyRm7Pn1XrJ9O5mPP8FPp/Rg3U1/JW/+fB+jNDWhNvOCb1OPiUhjYDbQDjjWG6qAiJwEfAa8CNyJe0DtYWAgkAi0VdXyI9zZXRQvAdaq6kkh2vd6s/Z2qDF1w9L1OfR//HMAkhKEubedSuumwW/8VF1VEmy8jdkVkdOAj6uw66eq2rOCcyQC83G5eb+gNsu5pspK8vLI/eADtrzwIgUrV5Zrb9K3L61uuJ4GnctNzmTqoNrMub4UuyKSCrwPHAn0UNXFQe3XAvcDpaPUP8HN4jAEaKyqhZWc+2ngclVtEKLNHpYwJo6c9+xcvlrl3ua8sfdB/LXP72rkOvH6gJqIpAH7V2HXHaq6ppLzjAQeANqp6oaA7ZZzTbVpURG5773H1kmT2Pnd9+XaG/U4hZZXXUVqt25IQiyMxjTRVqdnYwDwhidMBU4B+qhqyPcmRKQB0BnIVdW1IvIB0EhVT9nL+Z8GLlPVcl08lniNiS/vfr+e6yYtAqBVkwZ8ceuppCRF/5dfvBa70RJQ7O6jqhsDtlvONRHZ+eOPbBrzAHlz55ZrS0xPp8UVl5MxdCgJDcr1b5k6rE4XuyKSALwG/BE4S1VnVPG4w3Fz7w5T1Vcr2a8pbhjDqlBFsSVeY+JLYXEJJz8wi425OwF47MKuDOjafi9HVZ8VuxUTkSRgAdBSVTsEtVnONRFTVfIXLWLLiy+y/ZPyZUNS69a0uPwyMoYOtZ7eOFHXi91ngGuAfwLvBjWvU9V13py6w4G5uBkYjgFGAR+q6gUB5xoBHAzMAtYDHYDSbb1V9bMQ17fEa0yceWLGTzz88Y8AdNs/nbeuPTHq17Bi1xGRi4ABuGFoa4E2wF+Ak4CLVPW1oP0t55qoyl+ylMynnmL77NkQ9P8qsWVLWl51Jc0GDCCxWTN/AjRRUdeL3VW4ojSUu1V1tIi0ASYCXYEmwM/AC8BjqloUcK4/ArfhittmQC7wBXCvqn5ZwfUt8RoTZzK37+KE+2dSUFwCwDvXncgR+6ZH9RpW7Doi0h24Dzc7TnMgD1gIPKiq00PsbznX1Iji7GyyJk9h6/jxFG/dWqZNGjQg/bzzaDn8GpJatPApQhOJOl3s+s0SrzHx6W+Tv+XNRb8CcM5R7Rl7fteont+K3fBYzjU1rXh7HpvHPkzWq6+V6+lNSEuj+eWX0eLSS0lo1MinCE04rNiNgCVeY+LTd2uzGfDUFwCkJCYwd9SptGwcvQdWrNgNj+VcU1uKt20j5+13yH7jDXb98EOZtsQWLWh57XAyzjsPSUnxKUJTHVbsRsASrzHxa+BTX/DtWreS2oi+v+O6Uw+K2rmt2A2P5VxT21SV7bNms2nswxSs+LlMW1KbNjTp04f0c8+h4aGH+hShqQordiNgideY+DV10a/cNPlbANo2bchnt/YiOTE6T2ZbsRsey7nGL1pcTM7b77D58ccp2rixXHvDLl1IP28wTfv3J7Fp0xBnMH6yYjcClniNiV8FRSWcMGYmmdt3AfDUn46i/xH7ROXcVuyGx3Ku8VvJzp1kTZzElueeozg7u1y7NGhA0/79aXXjjSS3ae1DhCYUK3YjYInXmPg29uMfeXzGTwAc27E5U645PirntWI3PJZzTazQwkJ2fPUV2W9NZdv06WhBQZl2SU4mY9hQWlx+OUnNm/sUpSllxW4ELPEaE99+y93JiWNmUlTifsbfv+FkurSL/C1KK3bDYznXxKLinBxypr1L9n//y65ly8q0SVoaGRdeSIs/X0pSq1Y+RWis2I2AJV5j4t/1ry5i2nfrAbjgmP14YPAREZ/Tit3wWM41sUxV2fbJJ2y8a3T5uXpTUkgffC7Nhw0jpWNHfwKsx6zYjYAlXmPi39ert3LuM/MAaJCUwPxRvcloFNl0Q1bshsdyrqkLtKiInHemsfWll9j100/l2pv2O5OmZ51F4169ducCU7Os2I2AJV5j4p+q8scnP2fJr7kA3HbmIVzT48CIzmnFbngs55q6REtK2D5rFpnPjmPn4sXl2ht0OZQWf76Mpqf3tfl6a5gVuxGwxGtM/fD6wrXc8t/vAWifnsqckb1ITAg/Z1qxGx7LuaYuUlXy5s4l86mnyf/mm3LtSa1b03rEzTQ96ywkITrTG5qyrNiNgCVeY+qHnYXFnDBmJlvz3BPX44YezemHtQ37fFbshsdyrqnrdnz1FZsee4yd332PFhaWaUs58ECaDxtGswFnk9CwoU8RxicrdiNgideY+uNfHy7j6dluBaXjO7Xg1au6h30uK3bDYznXxIuirVvJmjCRrIkTKc7JKdOW2LIlLa++mvQLzifBhjdEhRW7EbDEa0z9sT47n5P/NYtibxqy6TedwsFtm4R1Lit2w2M518Sb4u15bHnuObImTKAkL69MW9I++9By+DWkDxqEJCf7FGF8iHbOtcEmxpi41C49lb5d2uz+evy8Vb7FYoyJD4mNG9H6rzfRedZMWt5wPYktWuxuK9qwgY133sXP/fqT8847aFGRj5GaQNaza4yJW/N/2cKF/54PQGpyIvNH9aZZWvV7XKxnNzyWc028K9m1i+zJk8kc92+Kt2wp05bSoQPNL7uMZoMG2vCGarJhDBGwxGtM/aKqnPnYZyzbuA2AO/ofyhUnd6r2eazYDY/lXFNflOTlsXXiJLa88AIlQWN6k9u3p9WNN9jsDdVgxW4ELPEaU/+8+uUaRr3p5szcv3kas0b0rPY0ZFbshsdyrqlvinNzyRw3juzJUyjZvr1MW1K7fWh13fU0O6u/zdO7F1bsRsASrzH1T35BMd3vn0FOvps26IVLjqH3oW32clRZVuyGx3Kuqa+Kc3LIfOZZtk6cCEFTlgG0e/BBmvbvZz29FbAH1IwxphpSUxK54A/77f765bmr/AvGGFMvJDZrRpvbbuWgmTPIGDq0XPv6W25h5YABZL/xJsXbtvkQYf1iPbvGmLi3dusOTnlwFqU/+jNu7sGBrRpX+Xjr2Q2P5VxjnOLcXDKffoatr7wCxcVl2hKaNaPF5ZeTPmggSa1a+RRhbLFhDBGwxGtM/XXF+IV88sNvAFxyfAfuHvD7Kh9rxW54LOcaU1bBmjVkT5lC1qRXKdmxo1x7k9NPp/XNfyNl//19iC52WLEbAUu8xtRfn/+UyZAXFgDQKCWR+f/XmyYNqzYNmRW74bGca0xoxdnZZE2eQvYbb1C4Zk3ZxqQkWl51JekXXEhym9b+BOgzK3YjYInXmPpLVenzyBxWbHJPSI/+YxcuPfGAKh1rxW54LOcaUzktLCT7jTfZ8tKLFK5eU6497dhjaTn8GtL+8AckKcmHCP1hxW4ELPEaU7+9Mm8Vf397KQCdWjbik7/1IKEK05BZsRsey7nGVI0WFLBt5iy2PP88O5csKdee2Lw5rW66kfSBA+vFtGVW7EbAEq8x9VveriK63zeDbbvcMp7jLzuWHr/b+wMhVuyGx3KuMdWjRUXkTJ1K9ptvkf/NN+XaEzMyaHPH7TTr39+H6GqPTT1mjDFhatQgicHH7Lv76/E2DZkxJoZIUhLpgwfTcdJEOr3/Pk37nVmmvTgri/U3j2Dl+RewbeZM+0Oyiqxn1xhTr6zMzKPXQ7MBEIHZI3rSoUWjSo+xnt3wWM41JnIl+flsef4FMp96KmR7mzvuIP3cc0hITa3lyGqO9ewaY0wEDmjZiJ4Hu6ELqvCfeat9jsgYYyqWkJpKq+uvo/OcT2lyxhmQXHYWmd/uvZefTjqZ9Xfcwc7ly32KMrZZz64xpt6ZtXwTf37pKwCaNExi/qjeNGpQ8ZPO1rMbHsu5xkRf4aZNbLzzLrbPnl2+UYQmffuSfv55NDrhhN25q66xB9QiYInXGANQUqL0HvspKzPzALh34O8Z0r1DhftbsRsey7nG1JydP/zAlueeJ3/xYgrXri3XnrTPPrQcfg3pgwcjCXXrjXwrdiNgidcYU+rFz1dyz7v/A+B3bRoz/aZTKuwFsWI3PJZzjal5qkr+t9+S+eyz5H06p1y7NGhA+gXn0+KSS0hu396HCKvPit0IWOI1xpTK3VlI9/tmsKPArVM/6YrjOKFzy5D7WrEbHsu5xtSu/CVL2frSS+S+917I9iZ9TqPtPfeQlJFRy5FVjz2gZowxUdC0YTLnHrVnGrKXbRoyY0wdl/r7w2j/8EP8bv48ml8yrNzDbNs+/oQVvU5l/R13ULDuV5+irH3Ws2uMqbdWbNrGaWPd234JAnNG9mLfjLRy+1nPbngs5xrjr6KsLPK+mEvWxInkL1pUpk2Sk2l2zjmkn3sOqUcc4VOEodkwhghY4jXGBBvy/AI+X5EJwNU9OjHqzEPL7WPFbngs5xoTO7KmTCFr0qvsWrasXFvT/v1pfcsIktu29SGy8qzYjYAlXmNMsI//9xtX/mchAOlpycwf1ZuGyYll9rFiNzyWc42JLapK/sKFbBr7SPme3gYNaNK3L63+ci0pHTv6E2BpLFbshs8SrzEmWHGJ0uPBWazLygfggXMP54I/7F9mHyt2w2M515jYpKrsWLCArNcms+3DD8u0ScOGZFx0Ea1uvIGEhg19ic8eUDPGmChKTBCGHb9njt2X56624swYE9dEhEbdu7Pvo4+w/3/G0+CQQ3a36c6dbH3pJZZ37cbGf9xLyY4dPkYaHdaza4yp97J3FND9/hnsLCwBYMrVx3PsAc13t1vPbngs5xpTN6gq22fP5rcxYyhcvaZMW9I++5Bx0UWkn3sOSS1a1Eo81rNrjDFRlp6WwqBueyZbH2/TkBlj6hERoUmvXhz43nu0vmUEkpKyu61owwY2jx3LTyeexK8jbqFgzZpKzhSbrGfXGGOAZRtzOePRzwA3tOHzW3uxT7NUwHp2w2U515i6qXh7HlvGjSP79dcpzs4u25iURPNLhtHquutISE2tkevbA2oRsMRrjKnMBePmsWDlVgCu69WZEacfDFixGy7LucbUbcXbtpH77rtkvTaZXcuXl2lLatWKZuecQ8urryIhrfz85JGwYjcClniNMZX5YPEGhk/8BoAWjVL44rZTaZicaMVumCznGhM/dixaxPqRt1K4dm3ZhoQEWo8YQfNhQ5GkpKhcy8bsGmNMDenTpQ3tmrmpdrbkFfDe9xt8jsgYY2JDWrdudHrvXdrcfjuJrVruaSgpYdO//sUvAweS8/bbaGGhf0FWwIpdY4zxJCUmcHH3PdOQjZ+3ynoljTHGk5CSQvOhQzjgjTdoNmhQmbaCFT+z/tbb+PnMfmRNmYIWFPgUZXk2jMEYYwJszXPTkBUUuWnI3rz2BI7u4KYhs2EM1WM515j4VpKfz+YnniRrwoRyxa2kptLiyitofvHFJDZrVq3z2pjdCFjiNcZUxYjXv+O/X68DYEDXdjx+0VGAFbvVZTnXmPqhKCuLrFdfJWv8fyjOySnTltCkiZun9/zzSdm3fQVnKMuK3QhY4jXGVMWSX3M464nPAUhOFFbc1x+wYre6LOcaU78UZWWR+cST5H7wAcVZWWUbk5Jo2rcP6edfQNpxx+4uaEOxYjcClniNMVV17jNz+Xq1S9arHzgLsGK3uiznGlM/aWEhuR98wKZHHqVoQ/kHfRsc1JkWV1xB0379kOTkcu1W7EbAEq8xpqrGfPADz376C2DFbrgs5xpTv2lBAdvnzGHry+PZsXBhufbkdu1oNnAg6RdeQHLr1ru3W7EbAUu8xpiqWL0ljzMe/Yz8wmL3tRW7YbGca4wptXPZMrKnvE7O1KmU7NhRtjEpiSa9etH0rLNo0uc0EhITASt2wxLLiTfgrxifIykrVuMCiy0csRoXxFZsd0xdzGtfrqWoxMVixW54YjnnRiKW/q9GW7zem91X7CjOySFr0iS2/ueV8uN6PV2WLwNsUQljjKkxUxet313oGmOMiZ7EZs1oOXw4nWd8QqubbiKpbdsav6YVu8YYEyRvV5HfIRhjTFxLSEuj5TVXc9DsWXR6/z2aDRxYc9eqsTOHICKDReQNEVktIvkislxE7heRJkH7dRWRD0Vku4jkisg7ItI5xPkaisiDIrLBO988ETml9u7IGBOPGjWIzvrudYmI/E1Epnn5VEVkdCX7DhSRRSKy08vnd4hIYi2Ga4yJIw06daLdmPvp9P57ZAwZEvXz13bP7gigGPg/4AzgGWA48LGIJACIyEHAZ0Az4GLgz0BHYI6ItA463wvAlcCdwFnABmC6iHSt8TsxxsStgd3akZRQ74bnXgm0BqZWtpOInA68AXwFnAk8BtwB3FfTARpj4luDTp1oe8ftUT9vrT6gJiKtVHVz0LZhwHigt6rOFJHngcFAR1XN9vbZF1gBPK6qI71tRwLfApep6kvetiRgKbBcVc8Ocf2YfVgiVgeZx2pcYLGFI1bjgtiKrT7OxiAiCapa4uXRQuBuVR0dYr9FQK6q9gjYdieu4N1fVTcGbI/ZnBuJWPq/Gm3xem92X3VLtKceq9We3eBC1/OV97l0DbnuwLzSQtc7bh2wBBgUcNzZuIQ8OWC/IuA14HQRaRDF0I0x9UiHFo14eshRpCYn1pseXlUt2ds+IrIf0BWYENT0CpCM6+k1xpiYEgsD00p7B37wPhcDBSH22wUcKCINVXUncBiwUlWDJmtjKZACdPb+XU5lS9T5LVZji9W4wGILR6zGBbEdm+Ew7/OSwI2qulJEdgBdQh0Ur9/TeL0viN97s/uqn3ydjUFE2gP3AJ+oaunSGsuBo0UkOWC/JrgkK0CGt7k5EGqCtq0B7cYYY6KnNK+Gyr1ZWN41xsQg33p2RaQx8DZQhHsIrdTjwHnAs944sCTgYaCx177Xt9oqEs/j7YwxppSInAZ8XIVdP1XVnjUVh+VcY0ws8KXYFZFUYBrQCejhjckFQFU/F5G/APcDl3mbP8E9xDaEPT23WUCHEKcv7VnYGqLNGGPqg7nAoVXYL3gY2N6U9uhmhGjLwPKuMSYG1Xqx6w1P+C9wDNBHVRcH76OqT4vIC7hxt7mqulZEPgAWqGqht9tSYJCIpAWN2+2CG/O7okZvxBhjYpSXE5fVwKlLn4M4DJhXulFEOgJpwP9q4JrGGBOR2l5UIgGYCJwKDFTV+RXtq6q7VHWpV+geDpyGm5e31DTc07/nBZw/CbgA+EhVd9XEPRhjTH2lqmuA73BzoAcagpsd54NaD8oYY/aith9QewpXnD4M5IlI94CPfcHNqSsi/xSR/iJymojcBnwOvKmqr5aeSFUX4aYde1RErhCR3rhpxw4A7gq8qHfOJ7wV1nZ4qwN1rI0brkxVV5TzKbbTRWSmiGwUkV0isk5EpohIyKet/eSttqcicq/PcfT04gj+yN770TVPRPqJyJyAlQkXisipPsc0u4LXTEXkQz9j8+I7UUQ+EpFNIrJNRL4Rkcv2fmTdIyLHiMhg4BxvUxcvRw0WkbSAXf8P6CEiE0TkUxHJB/4BrMLNhFOVa8Xs6pcisp+I/FdEcryfkzdFZP8qHHeMiPxbRJZ5v2fWiMhEETmgNuLem3DvK8R5bvN+Pj+viTirK9L7EpFDReR1EckM+D18Y03GXMW4wr4vEdlfRMZ7/wfzReRHEblXRBrVdNxViC3sekxEEkRklIisErd643cicm6VLqyqtfaBS4Zawcdob582uDG6mbjpxv4H3AwkhThfKjAW2AjsBBYAPUPs1xP4DXgfmO5dr2Nt3nsFr8d8YAqul6QHcBOQ7W1P8Dm2i4AHcQt89ACG4t7CzAU6+P3aBcW5wfue3utzLD29OK7HzRdd+nFMDLxOV+N63h4B+gCnA7cCZ/kcV5eg16o78FfvdbzW59iOAPKBWcAA73Ub58U23O/vaQ3c78uV5OeOQfte6OXnEmAzMAlYDPzKXUdLAAAPuElEQVQMNKrCtSZ6ue5KoDfwpvdad/X5NUgDfsJNrTbQ+75X6b6Ah4AvgGu9nPkn3JSaW4D96up9BZ2nE7Dd+336eQz8n43ovnDDKXOBd7zjewFXAX+rq/cFNAJ+BFYCl3j3NNL7+ZocA9+znoRZjwH/9PLOCO++xnk5qN9ej/X7xmvpxU0I+PcV1XlxaziuViG2DfPiO9Xv+ELEdrAX281+x+LFk4H7Q+ciYqvYPc3v1yYoro5eorvJ71iqGO8LXkJr7nMc9+HG/zcO2j4Pt/CN76+Vj6/Njbg50TsHbDsAN7tOpYUCcKT3c/LngG1JuGkn36nD9xUqn3fwfhnfU1fvK+g8070CYzaxUexG8v1KwHWmveX3fUT5vvp6P199g7aP8Y5P8/newqrHcEuZ78Kt7Bi4fQbw/d6O93We3dqiVVgZyA9atRXlYskW73ORr1Hs8QCwRAOGt5iQLsP9wn3W70D2xnu7/Dxgmqr6/WR/Cq43PD9oew4+z1EeA84G5qvq7geBVXUlrmdzQBWOjdXVL8O+r1D5XFVX43q+/c7nkXy/ABCRPwFHAaNqJMLwRHJfPXEzloytsejCF8l9lQ4lyg3ano3LW75OBxhBPXY67t6CV2+cABy+t+FC9T1hx6LgFeV8JSKJIpIiIgfh/qLfCPheXIrISbhe8L/4HUsIE0WkWES2iMikcMbFRdlJuCfzLxSRn0WkSERWiJviL9YMAprgphr028ve58dFpJ2IpItI6dvuj/gXVkw4jKBV1DxLqWAVtaBjV2rlq1/6JZL7KkdEDsX1SPmdzyO6LxHJwP2fHxkDf4QGiuS+TvI+NxSR+SJS6I3Nf1zc9Kh+iuS+PsENgXhARLqISGPv2YwbgWdVNS+6odaaw3A9u8EzbZXOEFPp6xILywUbj4ReUc5vC4CjvX+vwA2v2ORjPIhICq7wfkhVl/sZS5Ac3MOXn+L+qu6Ge5hnnoh08/F1a+d9POjF8zOu9/RJEUlS1cd8iiuUYcAmYuCpflVdIiI9gbdw4zDB9Uheo6qv+RZYbKhsBctQc/BW9djSdr9Ecl9liJsd6Flcz+4LkYcWkUjv60HcONCXoxhTNERyX+28z5OBJ4HbcGN47wH2w/3h7Zew70tVd3qdQW+wpxAEeB64LmoR1r7mQLZ6YxcCVClvWLEbI6TiFeX8NhRoinswYQTwsYicpKqrfIxpJO7hxH/6GEM56mYIWRSw6VMRmQN8CdwA3OFLYO4dnCbApar6prdtpvcE7CgReTxEAql1ItION8XgY97b2n7HcxB7fmFcgxvOMAC3uuNOVZ3oZ3wm5j0JnAD0V9VQhUudICIn4/4IPSoW8kQUlb6zPUFV7/T+PVtEEoExInKoqvrdI19tItIQV8C3xv3+XgMcC9yJqy+G+xedf6zYjQFSyYpyfgv4YV8gbmGPVbi/gK/xIx5vSMDtuIHtDYLG9zUQkXRgm6oW+xFfMFX9RkR+BP7gYxhbgIMov3zsR8AZwD7A+toOKoQhuF9AsTCEAdwDaoW4GStKF7OZISItgMdE5NVYfR6gFmQRuoepoh6p4GNjdfXLSO5rNxEZg3uq/xJV/ShKsUUikvsah+uZXuflV3C1Q6L3db76N699JPdV+gxKqLw4BvfOnF/FbiT3dTluPHJnVf3Z2zZHRHKAf4vIs6r6XdQirT1ZQLqISNAfXVXKGzZm12dSdkW5fhpiRblYoarZuKEMfo6p6wQ0xA1Kzwr4ANfznAUc7k9olfKzR2TpXtpjpWC7BPguhhLx4bh4CoO2fwm0wPWc1FdLcWPognVh76uoLQUOkLJz95Ye6/fql5HcFwAicjtuWr8bVPWVKMYWiUju61Bc50Zgvj0RN01gFv72FEb6/7AyfubFSO7rcCAroNAt9aX3uSrLiMeipUAD4MCg7aVjdSt9XazY9ZFUY0W5WCAibYBDcGM+/fItbn694A9wBXAvYmipaBE5Bjdl25d727cGveV9Pj1o+xnAOlXdWMvxlOO9Tl2InV5dcA9jdvXGiAc6Djevdyw9qFPb3gG6i0in0g3esJgTvbbKxPLql5HcFyJyA3AvcLuqPllDMYYjkvsKlW+/wz1A1QvXWeOXSO7rA9wDT6HyIoCfz81Ecl8bgQwRCe6UOs77/GuUYqxtH+LeaQu1euMSb7aKivk531ptfuAWRxiMW3JYcX+NDsYNG/ArptJY7qX8xPr7+vx6vQX8HTdGsRduUYJluOlLfuf39zNEvLEwz+5E73t5Du4PmJtxi6OsAVr6GJcAM3Fv212Dm4fxOe81u9Tv750X4+NeImvtdywBMQ32XqPp3s9BX9w4TAXG+h2fz69NI9wflYu91+ZsXAH0CwHzEuOGKxQBdwYd/xquV/AK3OwW/8X9AXFUXb0v3EIbJbgiKjifd6mr91XB+WYTG/PsRvr/8C5v+3245wVuw43Nf7mu3hduXvVc3AOFpYtK3OJtW4jPC1Z5Me61HvPu64Wg48Z4eeJvuKEaz3g/c3tdHMnXG67lF7eilYFm+xjTqkriGu3z63Ur8DWuuN2Bm/B9HDGwGEcl31+/i91RwPe4WRkKgbXAv4F9YuD1aYpbrvs33FvF3wN/8jsuL7Zk3BPr0/yOJURsZ3q/2DcD23DvLFwLJPodm98fwP64B/hyvddmanB+8H7xlstnVHH1y7p0X1S+At3sunpfFZxrNjFQ7Ebh/6HgCqcVXl5cjZuNIbmO31cX3Oqsa3HF+4+4Ff4y/L4vL769/px4X78cdFwi7kHv1bhe+e+BwVW5pngnMMYYY4wxJu7YmF1jjDHGGBO3rNg1xhhjjDFxy4pdY4wxxhgTt6zYNcYYY4wxccuKXWOMMcYYE7es2DXGGGOMMXHLil1jKiEiHUVEReTSKJ7zZRFZFa3zGWNMbRCRgSLyt2oec7SI7BCR9jUUz28i0jja5zbxxYpdYyq3ATgeeM/vQIwxxmcDcYswVMeDwIuqWhPL1L6Ny9G31MC5TRyxYtfEDBFp4HcMwVR1l6rOV9XNfsdijDF1iYgcjVuu9pmaOL+6VbH+DVwnIg1r4homPlixa6JORC4SkWUislNEFovI2SIyW0RmB+zT0xsecI6IPCcim3FL2Za2nyEi80QkX0RyRGSqiBwcdJ1VIvJyiOuriIwO+Hq0t+1wEZnlvaW2QUTuEZFKfwZCDWPwhiGsE5FuIvKZd76fROSaEMf3FpFvvNfiZxG5uoLrpInIAyKyUkQKvM+3B8YnIq+ISLaIdAjY1k5ENovI65XdhzHGRMLLtZcA7b2cqFUYjnUF8L2qLg061yoRmSAiF4rIDyKSJyILReSkoP3+ICIfi8gW73fBLyLydNA1pgDpwDkR3aCJa1bsmqgSkT7ARGAZLvk8BDwK/K6CQ57ArU8+FLjUO8cZuGED24ELgOHA74HPIxz3NRX4BPdW3CTg78CdYZ6rqXeOCcAA4CvgGRHpVbqDiBwKvI9bm/xC4P+Am4DegScSkSRgOu4Xw2PAmcDzXnwPBux6LbAFmCgiiV4h/AqwA7gyzPswxpiq+Acun23GDe06Hhi0l2POAD6roO1k4GZcnrsASATeFZF0AG8c7nSgGPe74UzgHiAp8CSqmgn84F3LmJCS9r6LMdVyN/A/YJD3FhMisgRYCPwYYv8vVfWKoG33Ar8AZ6pqkXeOed7xN1P9MWOlnlPVMd6/PxKRpsDNIvKoqmZX81xNgGtVdZYX3xzgdOAiYJa3zx3ANqCvquZ5+80FfgbWB5zrIuAkoIeqzvG2zRARgLtE5AFV3aSq20TkIuALXJG+C+gB9AwjfmOMqTJV/dl7B65AVefvbX8RaQN0BL6rYJemQFdVzfL234jrNOiH60g4BMgARqrq9wHHvRziXIuA7lW7E1MfWc+uiRoRSQSOAd4oLXQBVPVrYGUFh70VdI5GwFHA5NJC1zvHSlyR1yOCEKcEff0a0BjXa1xdO0oLXXBje3HF+P4B+xwPvF9a6Hr7rcXdR6AzgNXAXBFJKv0APgKSCUjiqvolrifkdtwfFv9U1c/DiN8YY2pSO+9zRc87zCstdD2Lvc+lOfQnIBsYJyJDRGS/Sq61OeB6xpRjxa6Jppa44mxTiLbfQmwD9yRtoAzcsIbg7QAbgeZhR1c+htKvwxkakRVi2y4g8CGJfUJcM1QcrYEOQGHQx5dee4ug/ScB6n08Va2ojTGmdpTmwl0VtG8N/MLrMNh9nKrm4B5uWw88DawRkSUicm6Ic+VTNvcaU4YNYzDRlIkr0lqHaGsDrAmxXYO+zvK2tQ2xb1vKJsidQErgDiISXBgGx/BL0NcANTElDriCvU2I7cHbtuB6vs+v4DyrSv/hjdMdD6zDDaV4FnswwxgTe7Z4nzPCPYGqfguc673TdQwwCpgiIkeq6pKAXZsHXM+Ycqxn10SNqhbjxuaeK96AU9g9/cwBVTxHHvA1cJ43LKL0HB2AE4DZAbuvpvwQhP6VnD64mLwQ9xDc4hD7RsM8oJ83NAMA7624E4P2+xDYD9iuqgtDfGQG7DsKN773T8BlwKCKZngwxpgo2wWkVnHfVbgOiU6RXlRVi7xxwn/H1S2HBu1yALA80uuY+GU9uyba7sKNNX1LRP6NG9owGjcEoaSK5/g7bjaGd71pZhrjxqfmAA8H7Pca8KKIPAK8CxyJN6NDBa70eka/wj1MdgUw2nu7rCbcC5yHexjuQVwv9GjKD2OYCPwZ91Daw7gHOlKAA4GzgYGqukNEjvOOv1tV5wF4r89YEZmjqj/U0H0YYwy4h4+bi8hwXMfGTlUN2VmgqgUisgA4NpwLichZwFW4WXRWAo2AG3AP/c4L2E+8awRPSWbMbtaza6JKVT8GLsb95f0WcCtuBoWNuGK1Kuf4ENdDm457qOxZ3NQyJ6lq4CwG43HF9TnANFwBW9lUOAOAPsA7wBBcMfqPKt5atXnFZz8gDZgMjMFNLTYjaL9CXOzP4ZL7+7gC+BJgLlDgzRwxyfv6voDDb8YNzZgkMbgohzEmrjyP62S4D/dMwbS97D8ZODXw3a1q+Ak3FvfvwAfAS0AR0EdV1wXsdwJuqMRrYVzD1BMS8NC8MTVCRPYFVuBmDqix4rKS64/GFcXJgTM8GGOMqTneH+nrcNM0TqihazwD/F5VT66J85v4YMMYTFSJSCowFrd4QyZuvNZI3MIHz/sYmjHGmFqkqrki8gAwUkQmapR710SkLe4dMFtQwlTKil0TbcW4WROexE2ZlYdbQec8VQ01nZgxxpj4NRa3Oto+lF1MJxo6AjcHLMZjTEg2jMEYY4wxxsQte0DNGGOMMcbELSt2jTHGGGNM3LJi1xhjjDHGxC0rdo0xxhhjTNyyYtcYY4wxxsSt/wcdX7T7TRPG7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(10,5))\n",
    "subplot(1,2,1)\n",
    "set_fig_properties([gca()])\n",
    "group_idx = range(1,9)\n",
    "plot(group_idx, temp_ave,linewidth=3,marker='o',markersize=10)\n",
    "xlim([1, 8])\n",
    "gca().set_xticks(group_idx)\n",
    "ylim([290, 310])\n",
    "gca().set_yticks(range(290,311,5))\n",
    "xlabel('group index')\n",
    "ylabel('T (K)')\n",
    "title('(a)')\n",
    "\n",
    "subplot(1,2,2)\n",
    "set_fig_properties([gca()])\n",
    "plot(t, Ein/1000, 'C3', linewidth=3)\n",
    "plot(t, Eout/1000, 'C0', linewidth=3, linestyle='--' )\n",
    "xlim([0, 1])\n",
    "gca().set_xticks(linspace(0,1,6))\n",
    "ylim([-10, 10])\n",
    "gca().set_yticks(range(-10,11,5))\n",
    "xlabel('t (ns)')\n",
    "ylabel('Heat (keV)')\n",
    "title('(b)')\n",
    "tight_layout()\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(a)** Temperature profile in the NEMD simulation. **(b)** Energies accumulated in the thermostats."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The figure above shows the temperature profile and the energies accumulated in the thermostats. The energy of the thermostat coupling to the heat source region is decreasing, because energy is transferred from the thermostat to the atoms in the source region. The energy of the thermostat coupling to the heat sink region is increasing, because energy is transferred from the atoms in the sink region to the thermostat. The absolute values of the slopes of the lines in **(b)** should be the same; otherwise it means energy is not conserved. The absolute slope is the energy transfer rate $Q=dE/dt$.\n",
    "- The thermal conductance in this system can be calculated as\n",
    "$$\n",
    "G = \\frac{Q/S}{\\Delta T}\n",
    "$$\n",
    "where $S$ is the cross-sectional area and $\\Delta T$ is the temperature difference, which is 19 K here (slightly smaller than the target value of 20 K). The calculated thermal conductance is about 10 GW m<sup>-2</sup> K<sup>-1</sup>. This is the classical value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.04638256513016"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "deltaT = temp_ave[0] - temp_ave[-1]  # [K]\n",
    "deltaT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.635182"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q1 = (Ein[int(ndata/2)] - Ein[-1])/(ndata/2)/dt/Ns\n",
    "Q2 = (Eout[-1] - Eout[int(ndata/2)])/(ndata/2)/dt/Ns\n",
    "Q = mean([Q1, Q1])  # [eV/ps]\n",
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.82366678786541"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l = gnr.cell.lengths()\n",
    "A = l[0]*l[2]/100  # [nm2]\n",
    "G = 160*Q/deltaT/A  # [GW/m2/K]\n",
    "G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Spectral Heat Current Results\n",
    " - The [shc.out](https://gpumd.zheyongfan.org/index.php/The_shc.out_output_file) output file is loaded and processed to create the following figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['t', 'Ki', 'Ko', 'nu', 'jwi', 'jwo'])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shc = load_shc(250, 1000)['run0']\n",
    "shc.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ly = split[5]-split[4]\n",
    "Lx, Lz = l[0], l[2]\n",
    "V = Lx*Ly*Lz\n",
    "Gc = 1.6e4*(shc['jwi']+shc['jwo'])/V/deltaT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArsAAAFTCAYAAAA5nMTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5xU9fX/8dd7O0vvvYMgiL1g+0aNsddoYolEU4z6Tf2ZHv0mtnSjKSYxGhMTY2KLiiVRrNhARQVkKdJ7WdqylIUt5/fHvbMMy+zu7O70Pc/HYx4Ld+7cOSt4OfuZ8zlHZoZzzjnnnHO5KC/dATjnnHPOOZcsnuw655xzzrmc5cmuc84555zLWZ7sOuecc865nOXJrnPOOeecy1me7DrnnHPOuZzlya5zzjnnnMtZnuy6rCcpX9JcSU+08vVdJG2W9ItEx+accy72fVrSq5LibvYv6UZJlZL6JCdKl6s82XW54IvAWODm1rzYzLYBdwJflTQkkYE555wD2nifDv0WqAZuSkRArv3wZNdlNUkFwP8BL5rZzDZc6i6C/x++n5DAnHPOAYm7T4cLE38GrpY0KFHxudznya7LducAA4EH23IRM9sC/Af4jKROiQjMOecckKD7dOhBoAD4QgKu5doJT3ZdtrsKqAP2qdeVdISk30sqk7RN0g5J70v6X0lq5FqPAZ2Bi5MbsnPOtStXEeM+HSGpg6RfSVolqUrSLEmTYp1rZrOAReE1nYuLJ7sua0nKAz4GzA0/3op2NXA+MAu4G/g70BX4PUF9bizTwq+nJD5a55xrf5q5T0c8ClwEPEJQptAf+Luk/9fI+dOAYZJGJDpel5s82XXZ7ECgG/BejOd+Agwxs8vN7Dtmdh0wBnieYCPa0IYvMLMlwBbg+CTG7Jxz7UlT9+mI4cBBZna9mX0FOBTYAPxU0oAY588Iv/q92sXFk12XzSIbFNY3fMLMVphZXYNjNcA9BH/vT27kmuujruucc65tGr1PR/mxmW2P/MbM1gC/AYqBS2KcH7mW36tdXDzZddmsR/h1a8MnJBVL+rakGWFfRgv7Of47PKV/I9fcDBRJ6pyEeJ1zrr1p9D4d5fUYx94Ivx4S47nN4dderQ3KtS8F6Q7AuTaoCr+WxHjuceAsYD7wT6AcqAGGAVcSrBjE0gEwYFciA3XOuXaqqft0xIYYxyKrt11jPNch/LqztUG59sWTXZfNysOvPaIPSjqKINF9Djg7upxB0iUEyW5jegBbw5IH55xzbRPzPt1AH2Blg2N9w68VMc6PXKs8xnPO7cfLGFw2KyNYhR3d4PjI8OuzDet2aWJDg6SOBDVgHyYsQueca98au09HOzHGsRPCr7NiPDcm/Or3ahcXT3Zd1goHQcwBjmrw1Irw6z6JraSJwJeauOQRQD4wNVExOudce9bEfTraDdHDfCT1B74O7CZoR9bQMcAeYHoCQ3U5zJNdl+0mAz0kHR517G2C1jSXSnpF0i8kPQa8BjzTxLVODb8+mZxQnXOuXYp1n462FJgTDpb4HcFqbh/g+2a2OvrEMCmeCDxvZr63wsUlrcmupEGSfidpmqSd4Y75YTHOK5H0S0lrJe0Kz/+f1EfsMtCfCSbzXBE5YGa1BOMp/0bwcddXCPo4XgXc1cS1LgfeM7P3kxWsc+kgabCkxyRVhBMFH5c0JI7XDZU0WdLy8N67UdJUSWfFONcaeRyanO/KZZH97tMNfIqgU84lBJ++rQeuNLNYA4AuJNigdk8S4nQ5SmaWvjeXTgIeJmg2nQ+cBgw3s2UNznsQOBv4NrAE+DJwJnCsmc1MYcguA0l6FPgfYFhrf9KX9DHgVeAKM0vE/HbnMoKkUoKVst3AjQT1k7cBpcDBZrajideOB64n+H9jFdCFYDrh2cBFZvZ41LkG3A/8qcFlZpuZ75pv5xJxnw6v8wpB68jx4cKGc81Kd7KbF9lAJOmLwL00SHYlHQLMBD5vZn8NjxUQFL0vMLPzUh64yyiSDiD4+/BtM/t1K6/xEtAdOMLS+T+Fcwkm6evAHcAYM1sUHhsOLAS+Y2Z3tPB6BQQfO880s3OjjhvBcIAbExa8yxkJuk+fSFCOdqGZebmZi1tayxhi7JSP5TygmmAFOPK6GuAh4HRJjfVLde2EmX1EUKJQ1cypMUnqQnADvdoTXZeDzgOmRxJdADNbCrwJnN/Si4X33wqCvtXOxaWt9+lQd+Bbnui6lsqGPrvjgaUxPgYrA4qAUeGvgfrVBddOSfpjG15+k6SExeKyh5nl8h/8eIINQg2VEdRKNktSHsHiSC+CmsoDCHbLN3SdpG8DtQQ75X9kZvtNx/L7dPvWxvs0km5PVCwue7TlPp0N3Rh6AFtiHN8c9bxzzrnYmrqHdo/zGr8g+IRtLcHeiUvN7KUG5/wD+F+CriZfAnoCL4d7M5xzLm2yYWW3VfzTaBctsmLrfy9cNF/Jj9uvCUrH+gGfBf4p6WIzq2/lZ2aTos5/XdJkgv6qt7F3QMA+cun/x1y9x+Ti95WL3xPk/vfVFtmwsruF2KsPkRXdzTGec845F2jqHhprxXc/ZrbKzGaY2TNm9mmCEoUmP0o2s0rgWZoeJuCcc0mXDcluGTA8bJ8TbRzBBJVF+7/EOedcqIygbrehccDcVl5zBsF+iXjk1jKTcy7rZEOy+zRQSNRGirD1zSXAFDPbna7AnHMuCzwFTJQ0InIgHN5zfPhci4Sb1U4AFjdzXheC4S7vtPQ9nHMukdJesyvp4vCXR4Rfz5RUDpSb2VQz+0DSw8CvJRUS9He8jmAi1mdSH7FzzmWVewmmCE6WFBkqcSuwkqgBEJKGEiSwt5jZLeGxmwjKHd4E1hHU7H4BOJpg4mDktd8imFb4CrAGGAp8Kzzf79POubRKe7ILPNrg938Iv04FTgp//TngxwQbHboRTAM6w8e6Oudc08xsh6RTgDuBBwABLwHfMLPtUaeKYJJl9Cd+7wPfAC4FuhIkvLOAE83szajzFhCMcb0wPG8bQYL8BTPzlV3nXFqldYJaMkT6N+ba9+XaJq+4IyosonZ7XPtxXDsRtXvZ2zKkUC7ep3N9J3wufV+5+D1Bu/i+Wn2fzoSVXeeSbsQ3H6Wmztixu4aOxf7X3jmXWLmWYETk4veVi98T5O73lQjZsEHNuTarqQtuAmsrdqU5Euecc86lki9xuZxXW7f3p92aOv/J1znn2pNde2qZPHM1D769grUVu5g4oic3nH0g/bt2SHdoLkU82XU5b8eemr2/3l2bxkicc86lgpkxb20lk2et5uF3V7J1Z3X9c8/MXst7y7fw1FdOoHfn4jRG6VLFk12X83bsron5a+ecc7lhXUUV/3xnBe8v38LKLTtZV1HF7pq6+ucPGdyNzx03jPEDuvDtx2Yzc+VWfvDEh9z72SPTGLVLFU92Xc6LXs3d7smuc87llKkflfOVB9+nssH9vU/nYj5+YB8+feRgDhuyd2L2nyYdwcm3v8oLc9fz7rLNHDWsR6pDdinmya7LedGruZ7sOudc7li4vpJrH3iPXdW1nDymN5cfM5ThvUrp37VDo513+nYp4QsnDOd3Ly/ivteXerLbDng3BpfzvIzBOedyj5lxw5Nz2FVdywWHDuAvVx3FJ8b1ZVSfzs22mJw0cSiF+WLK3HWUV+5OUcQuXTzZdTlvx56oMoYqT3adcy4XTJm7nneWbqZnxyJuPu+g+uED8ejTpYT/Gd2bOoPnytYlMUqXCTzZdTlvnzKGPZ7sOudcLvjrm0sB+PLJo+haWtji1599cH8Anp29JqFxuczjya7Lefu2HvNk1znnst38dduYvmQzHYvy+dSRg1p1jVPH9aUoP4+3l25m8449CY7QZRJPdl3O2xnVjWHXnromznTOOZcNnvwgWI09/7CBdC5p+aouQJeSQo4c1h0zeGvxxkSG5zKMJ7su5+2p3Zvg1tZ5suucc9nMzPjPh2sBOO+QAW261gmjewHwxkJPdnOZJ7su59XU7h0RXO3jgp1zLquVrdnGis076dWpuM1tw04c1RuA1xduxMz/fchVnuy6nBe9mltT6yu7zjmXzV6ZvwGA08b3JT8v/g4MsYwb0IVupYWs3rqL5Zt2JiI8l4E82XU5L3o1t9ZXdp1zLqu9GdbXnjiqV5uvlZ8njh/Za5/rutzjya7LedGrudW1nuw651y22rWnlveXb0WCY0f2TMg1jxkRlELMWLYlIddzmceTXZfzaqJWc2t8g5pzzmWt95ZvYU9tHeMHdKFbaVFCrnnk0CDZfXfZ5oRcz2UeT3ZdzoveoFbjK7vOOZe1IqUGx41sewlDxJh+nelcXMCqLbtYV1GVsOu6zOHJrst5+67serLrnHPZ6q1FkWQ3MSUMENTtHja0OwAzlvvqbi7yZNflvOiaXe/G4Jxz2aliVzUfrq6gMF8cPbxtLccaOjKS7Hrdbk7yZNflPF/Zdc657Dd9ySbqDA4b3J3SooKEXvvIYb6ym8s82XU5b59k12t2nXMuK01bvAlIXBeGaIcO7kZ+npi7Zhvbd9ck/PouvTzZdTlvnzIG78bgnHNZ6c2wXvf4BPTXbai0qICDBnShzmDmiq0Jv75LL092Xc7zMgbnnMtuG7ZVsXDDdjoU5nPo4G5JeY8jvAVZzsqaZFfS8ZKmSNogqVLS+5I+n+64XObbd4OaJ7uu/ZE0WNJjkiokbZP0uKQhcbxuqKTJkpZL2iVpo6Spks6KcW6JpF9KWhueO03S/yTnO3LtzVthCcNRw3tQVJCc1OWosG73veW+SS3XZEWyK+lg4EWgELga+CTwLnCfpOvSGZvLfD5UwrVnkkqBl4GxwJXAJGA08Iqkjs28vBOwEbgROAv4AlAJPCvpkw3OvY/g/vxD4BxgLfC8pEMT9K24dqy+hCEJ9boRR4TJ7vsrtnjnnhyT2O2MyXMpkA+ca2bbw2MvhEnwZ4E/pi0yl/F8qIRr564GRgBjzGwRgKTZwELgGuCOxl5oZmUECW49Sc8CS4HPAY+Hxw4BLgc+b2Z/DY9NBcqAW4DzEvstufbEzJJarxvRp3MJQ3uWsnzTTsrWbOOQJJVLuNTLipVdoAioBnY1OF5B9nwPLk2iV3Or/ad11/6cB0yPJLoAZrYUeBM4v6UXM7Magntv9Jb18wju0Q83OO8h4HRJxa0L3TlYXL6dNRVV9OxYxLj+XZL6XhOHByvH05ZsSur7uNTKlkTx/vDrbyUNkNRN0tXAx4E7Y71AUqMP175ElzHU+ga1dsfvBYwH5sQ4XgaMi+cCkvIkFUjqJ+mHwAHAXQ3eY6mZ7YzxHkXAqEau297/bFwcpn4UrOqeOLoXeXnJ/bsRaWsWaXPmUiPZ94KsKGMwszmSTgKeAP43PFwNXGtmD6UtMJcVoksXqj3Zde1PDyDWjpvNQPc4r/EL4Jvhr7cDl5rZS3G+R+R551rl9YXlAJw4unfS3yuS7L67bDPVtXUU5mfLmqBrSlb8KUoaDfybYJXgXOBU4G7gbkmfifUaM2v04dqXah8X3K75vSAhfg0cRXD//S/wT0nntPWi/mfjmlNVXcv0sKTgxAOSV68b0bdLCSN6dWTnnlpmr6pI+vu5QLLvBVmxsgv8hGAl9xwzqw6PvSSpJ/AbSf8yM89iXEzRpQt1BnV1lvSPwpzLIFuIvYLb2GrsfsxsFbAq/O0zkl4FbgeeiXqPoY28B+xd4XWuRV5fuJGq6jrGD+hCn84lKXnPiSN7smTjDqYv2cQRQ+P98MNlsqxY2QUmALOiEt2Id4CeQJ/Uh+SyRcNBEj5YwrUzZQQ1tQ2NA+a28poz2LcOtwwYHrY5a/gee4BFONcKT89aA8A5Bw9I2Xse53W7OSdbkt11wKGSihocPwaowlcNXBMa9tb1TWqunXkKmChpROSApGHA8eFzLSIpDzgBWBx1+GmCPuifijqvALgEmGJmu1sTuGvfdu6p4cV56wE45+D+KXvfY0f0RIJ3lm5m++6a5l/gMl62lDHcBTwKPC3pDwQtyM4DLgPuNLM96QzOZbaGvXWr6+roQH6aonEu5e4FvgJMlnQjYMCtwErgT5GTJA0lSGBvMbNbwmM3EZQivEmw6NCPoO/u0QR9dQEwsw8kPQz8WlIhQR/e64DhQMx9Fa59MzP+Pm05D0xfTmlRPl84YTjnHTJgn933T3ywmp17ajl8SDcG92j4oUHy9OxUzBFDujNj+RZe+6icsyakLtF2yZEVK7tm9hjB9J5i4M8Em9VOAL4MfDuNobksUN0g2fXBEq49MbMdwCnAR8ADwIMEyegpUUN6AEQwvCf634X3gYOA3wFTCLoyVAEnxuiE8zngr8BtwLPAYOAMM3s/0d+Ty373vr6EHz1VxqIN25m9qoKvPzSTbz46i6rqWiDYmPanqUsAuOr44SmP7xPj+gIwpWxdyt/bJV62rOxiZv8l2AXsXIvUhmUM+Xmits58ZLBrd8xsBXBRM+csI0h4o489RZylDma2C7g+fDjXqDVbd/HL5xcA8OMLD6LO4CfPzuPx91ezuHwHv/rUwdz/1jJWbN7J6D6dOPOgfimP8bTx/fjpf+fz8vwN3oIsB2RNsutca0VWcksK8tixp9ZXdp1zLo3ueW0J1bXGuYcM4DPHBE08jhjSnav/PoNZK7dy6h2vAVCYL37yyQlpSTSH9+rI6D6dWLhhO28s3MjJY30ffDbzH1Vczot0XygpDOp0Pdl1zrn0qK6t48mZqwG49mP1eyYZN6ALT33leM47ZACdSwoY3acTf77yKI4alr55JBcePhCAh95dkbYYXGL4yq7LeZGyhUiyW+1lDM45lxZvLNrI1p3VHNC3E+MHdN3nuZ6divntZYelKbL9XXz4IH415SNemreBDZVVKevz6xLPV3Zdzous7BYXBH/d67z1mHPOpcUzs9YCqe2b21p9upRwytg+1NQZj85Y1fwLXMbyZNfltLo6IzJtsCA/2HtT66NInXMu5cyM1xaWA3BGGjadtcakiUFN8V/eWMquPbVpjsa1lie7LqdFEtv8PJEX9m/0oRLOOZd6i8t3UF65m16dihndp1O6w4nLiaN7ccigrmzasYcH316e7nBcK3my63JaJLHNl8jPC5JdL9l1zrnUm7Z4IwDHjuy5z/CITCaJr586GoC7py6msqo6zRG51vBk1+W0unBlNy+P+mTXyxiccy71pi3ZBATjeLPJyWP6cMTQ7mzcvoffvbwo3eG4VvBk1+W06JVdL2Nw2UDSREk3SXpO0mxJCyVNk3S/pM9J6p7uGJ1rjQ9WbAXg6OHpayfWGpK46dzxSEHt7uLy7c2/yGUUT3ZdTouULOTliYI8T3Zd5pJ0paQPgbeA/weUAguBt4EtwDEE49JXh4lv6meoOtdKGyqrWFtRRafiAkb06pjucFpswqCuXHLkYGrqjFufmZvucFwLebLrcto+G9Q82XUZStJs4GfAf4AjgG5m9j9mdpGZXWFmZ5nZgUAP4GqgDzBX0iXpi9q5+M1ZXQHAQQO71N+Ls823Th9D55ICXl1Qzotz16c7HNcCnuy6nLbPBrWwjKHOa3Zd5rkPGG5m3zWzD8xi/yU1swoze9DMzgImAltTGqVzrTR7VZDsHjyoW5ojab1enYq5/hMHAHDT02XeiiyLeLLrctreDWp7uzH4yq7LNGb2GzOrauFrZpnZ88mKyblE+jBMdicM7NrMmZlt0sShjOvfhVVbdnHXKwvTHY6Lkye7Lqfts0HNuzE451zKmRmzV0dWdrM72S3Iz+O2Cw8C4J7XlrBoQ2WT55sZz8xew01PlTF55mqf4Jkmnuy6nFaf7OaJcICa32xcRpP0F0nfauS5EZL+kuqYnGuLDZW7Ka/cTZeSAob0KE13OG12+JDuXHb0YKprjf97soxGqo6oqzO+8fBMvvLPD7j/rWV8/aGZ3OKb29LCk12X0yJlDBJexuCyxVXAzyU9LKm4wXO9gStTH5JzrbdwfdCqa0y/zlkzTKI53zl9LD06FjFtySYmz1wT85yf/ncek2euoVNxAVcdN4yigjzuf2sZbyzcmOJonSe7LqdF8tp8r9l12eVG4FTgVUm90x2Mc20R+ah/VJaMCI5H945FfO/MsQDc9uxcKnbtO1ntL28s5d7Xl1KQJ/406QhuOm88X/94MInt1y9+lPJ42ztPdl1OizUu2Gt2XRZ4CTiWoNXYO5LGpTke51ptcfkOAEb2zp1kF+DiwwdxZDhZ7av/+oDq2qCx+2PvreLWZ4NyhV9cfDDHj+oFwFXHDaNzcQEzlm9h3tptaYu7PfJk1+W06G4MPkHNZRMz+4hgkMQS4C1Jp6c5JOdaZdGGoIxhZA6t7ELw78qvPn0IPTsW8dpH5Vz0x7f4yj/f51uPzsIMvn/mWD55+KD68zsWF3D+YQMAeHpW7NIHlxye7LqcFmtl1/vsumxhZluB04FHgKeBL6U3IudablE4XndUjq3sAgzt2ZH7P3c0fbsUM3tVBc/MXkt+nrjhrAO55mMj9zv/rAn9AXiubF2qQ23XCtIdgHPJFEl28/L2DpUIP2lyLiuYWQ3wJUnzgF+mOx7nWqJiVzXllbspKcxjYLcO6Q4nKSYM6sqL13+M5+aso7KqhpPH9mF4IyORjx7Wg64dCllSvoOVm3cyOAe6U2QDT3ZdTqurHxdMfZ9dbz3mMtzJwLyGB83sTklvA6NTH5JzrbM4XNUd0atT1o4JjkfnkkI+deTgZs8ryM/jmOE9mDJ3PdOXbPJkN0W8jMHltOgyhgLfoOaygJlNNbOYnerN7C0z+1uqY3KutRbnaL1uW0wc0ROAaUs2pTmS9sNXdl1O22eDWpjs1vjKrsswkn7YgtPNzG5NWjDOJVAu1+u2ViTZfXvJZswsZ3oPZzJPdl1Oi9Tn5mtvza6XMbgMdFOMYwbE+lfQAE92XVZYvCFsO9Yndg1rezS2X2e6dihk9dZdrK2oYkCO1jJnkqwqY5B0lqTXJG2XtE3SDEmnpDsul7n22aDmQyVc5ips8OhAkOgeE+O5opZeXNJgSY9JqgjvnY9LGhLH646UdI+k+ZJ2Sloh6UFJw2Ocu0ySxXhc0NJ4Xe6I1Ozm0kCJtsrLEwcP6grA7FUVaY6mfciaZFfSNcBk4D3gQuBTwKOAV3e7RtVvUNPePrveesxlGjOrjX4ANeFTtQ2fC5+Pm6RS4GVgLMGo4UkEm9xekdTcctulwHjgt8CZwPeAw4EZkmLtxnmeYBhG9GNqS+J1uWN3TS0rNu8kTzCsp6/sRoskux+u3prmSNqHrChjkDQM+DXwbTP7ddRTz6clIJc16jeo5Yn8vH2POddOXA2MAMaY2SIASbOBhcA1wB1NvPbnZlYefUDSm8DS8LoNa403mtn0RAWeS8yMaYs3MXPVVvp1KeGMg/pRWpQV/wS32vJNO6mtM4b2LKWkMD/d4WSUCQO7Ab6ymyrZ8n/a54E64O50B+KyS22MDWrejcG1M+cB0yOJLoCZLQ2T1vNpItltmOiGx5ZLKgcGJiPYXLR15x6++cgsXpq/of7Y7c8v4M9XHsW4AV3SGFly1Xdi8M1p+4kuY/BNasmXLWUMJwDzgUslLZZUI2mRpC839gJJjT5c+1FX33qM+tZjvkGtffF7AeOBOTGOlwHjWnoxSQcCfYjRCxg4N6zt3S1penP1uu3hz2ZbVTVX3Pc2L83fQNcOhVx57FAO7N+FNRVVTLrvbdZW7Ep3iEkTGRPs9br769+1hF6diqjYVc2qLbn7dyBeyb4XZMvK7oDw8UvgB8BigprduyQVmNlv0hmcy1z7lDHIW4+5zCRpRINDkc98B0rar6jPzJa04PI9gC0xjm8GurfgOkgqIPiErRy4r8HTTwPvEpQ49AW+AjwhaZKZ/aMl75MrzIzrH57JnNXbGNqzlAe/eAyDupeyu6aWz9//Lm8u2sTX/vUBD3/p2JwcuBDZnDayt9frNiSJMf06s3HRJhasq/ThEkmWLSu7eUBn4Bozu9fMXjaz64DngO8rRupvZo0+XPtR32dX8glq7VSW3AsWEdTQRh7zw+NPNjgeeaTLXcBxwBVmtk8CbWZfNbO/m9nrZvYY8HFgBvDTxi6WJX82rfb3act5cd4GupQU8MDng0QXoLggn99eehi9Oxfz7rItPP7B6jRHmhyLvBNDk8b0DUpYFqyPOUOmXUn2vSBbVnY3EewefqHB8SnAGUB/YE2qg3KZr77PbtTKrtfsugz0eYL+ucmwhdgruI2t+MYk6WfAl4ArzWxKc+ebWa2kR4GfS+pvZmvjfa9csHB9JT/+T1Dp8bOLDmZIz31X7np2Kub7Z47l+kdm8Yvn5nPOwf1zahNXXZ3t7bHrNbsxje3XGYAF6zzZTbZsSXbLgIlNPF+XqkBcdom5Qc3/trgMY2b3J/HyZQR1uw2NA+bGcwFJNwDfBb5qZg+0IoZ29ROmmXHDE3PYU1PHp48cxFkT+sc878LDBnLfG0spW7ONh99dyZXHDUttoEm0dlsVu6pr6dWpiG6lLW4N3S6M8WQ3ZbKljOGJ8OvpDY6fAawys3Upjsdlib0b1PYOlfA+uy7TSFoi6ZAkXf4pYGJ0XXDYzvH48LnmYvsacBtwg5ndFe+bhvW9lwAr2ts9+rH3VvHOss306lTEDWc1vgdQEl89ZTQAf3x1MbtrWtRCOaNFNqeN8FXdRo3u2wkpqG3eU+OrMMmULcnuf4BXgD9JulbSaZLuBU4D/i+9oblMFmuDmvfZdRloGFCcpGvfCywDJks6X9J5BAN6VgJ/ipwkaWjY6eaHUccuJehx/hzwsqSJUY9xUeddJukhSZ+VdHL4ulcIBlB8N0nfV0basmMPPwnLF248exxdSwubPP+0cX0Z268z67ZV8fj7uVO7u9g7MTSrtKiAIT1Kqakzlmzcnu5wclpWJLsWVChfADwE3Aw8QzBG8zNJ/vjPZbnaqA1qPi7YtUdmtgM4BfgIeAB4kKBjwilmFv0vrAi6QET/u3BGePwMYFqDxx+izltK0I7slwR7Ke4GdgNnmNlDif+uMtfP/jufLTurOW5kT84/dECz5+flietOGgnAPa8tyZn706Jy77EbjzF9vZQhFbKlZhcz2wZ8OXw4F33h0CcAACAASURBVJf6MoY8PNl1mS5pfzHNbAVwUTPnLCNIbKOPXQVcFcf1pxMk1O3au8s28/CMlRTl53HrBQfF3SP07An9uX3KApZu3MGUsnWc2UiNbzbxld34jOnXmSlz19eXfbjkyJpk17nWiKzs5vsENZf5bpa0MY7zzMyuTHo0rkWqa+u44YkPAbj2pJEtWtEsyM/j6hNH8MPJZdw9dTFnHNQv6wdreI/d+ET+nkT+e7nk8GTX5bTIym6e9tbsep9dl6EOJfjovzn+FzgD3f/mMj5av52hPUv537AsoSU+dcRgfvPiQmatqmDa4k0cN6pXEqJMja0797Bx+x46FOYzoGuHdIeT0eqT3bBNm0uOrKjZda619tmglrfvMecyzAVmNjyOR8Npay7NNm7fzW9fCmZ93HTu+Fb1y+1QlM9VYeuxP05dnMjwUi6ySjmid8ecnAyXSCPCle+lG3f4v01J5Mmuy2m14b0jTyLPh0o455LgV1MWULm7hpPG9ObksX1afZ1Jxw6ltCif1xduZM7qigRGmFqRVUqv121ex+IC+nctYU9tHau27Ex3ODnLk12X0+r2Wdn1MgbnXGLNWV3BQ++upCBP3Hh24z1149GttIjLjh4CwN1ZvLpbPybYOzHExet2k8+TXZfTojeoRZLdGk92XWbyz3uz0I+fnYcZXHncsISsZH7xxOEU5ov/fLiWsjXZubob6cQw0ld24xIpZfC63eTxZNfltNq6/fvs+gQ1l6GelHSPpDMl+XzVLDBt8SamLdlE55ICvhZOQmur/l078JljhlJncPPTc7EsvF/Vr+x6shuXyMquD5ZIHk92XU7bp8+uT1BzmWsAwcCcwQTj0cslPRJOJuuS3tBcLGbGnS9+BMDVJ45odlJaS/y/Uw+gR8ci3lm6mcfeW5Ww66ZCVXUtKzfvJE8wtGdpusPJCt6RIfk82XU5LXqCWn2fXR9B7jKMma0zs7vN7EygN3ANUAv8kSDxnSLpOknNj+RyKTFt8SbeWbqZrh0K+dzxwxJ67a6lhfzgrAMB+NFTZSzJolrOZZt2UGcwtGdHigta3pWiPRrZJyxjyKI/52zjya7LaTH77Gbhx4Ku/TCzSjN7yMwuI0h8zwcWAzcCKyW9I+n7aQ3Scc/rSwD44gnD6VySuFXdiIsOH8g5B/dn555avvTAe2zesSfh75EMkdVJHyYRv35dSigtymfTjj1syZI/52zjya7LaZGKhegNal7G4LKFmVWb2XNmdp2ZDQSOB14GJqU5tHZtcfl2Xl1QTnFBHldMHJqU95DETz45gQP6dmLRhu1c9dd3siIRWuSb01pMktftJpknuy6nxRoX7Cu7LtNI6hHPeWY2HZhpZm3rceXa5G9vLQPgwsMG0r1j8vYSdikp5IEvHMPgHh2YvaqCi+5+i5WbM7sX694xwZ7stsRI78iQVJ7supwWXcZQEGk9VuvJrss4L8SzEU3SVcADyQ/HNaaqupYn3l8NBO3Gkq1vlxIeveY4xvbrzJLyHVz4h7cyeuBEZGXXOzG0zAjvtZtUnuy6nFYb1Y3BJ6i5DDYceE5SoxmCpC8B9wFTUhaV288Lc9dTubuGgwd15cD+qWmU0a9rCY9ceyzHj+rJxu27ueze6by3fHNK3rsl6uqs/mP4kb082W2J+l675b6ymwye7LqcFt2NwSeouQx2OjAOeFZSh4ZPSvoacDfwFHBBimNzUR5/P2gF9snDBqb0fbuUFPLXq47mrAn9qKyq4Yo/v8OslVtTGkNzVm3ZRVV1HX06Fye0FVt74DW7yRV3sitpmKRLJV0v6QZJ10g6SVJJMgN0ri32HRccHPOVXZdpzOxd4CzgMOApScWR5yR9G/g18ChwsZlVpydKt3H7bl5buJGCPHHuIanvAldUkMdvLz2MCw4dwK7qWr7wt3czqoZ3UXklAKP7+qpuSw3v1REJVmzaSbX3x0y4JpNdSd0kfUvSfILWN/8EbgduJej/+DKwNWx+flKyg3WupfbZoCZf2XWZy8zeAs4FjgMel1Qo6UfAz4EHgcvMrDadMbZ3L8xdT22dceLoXvTsVNz8C5KgID+PX37qEE4Y1YuN2/fwlX++nzHJ0cL1wark6D6d0xxJ9ikpzGdA1w7U1BkrMugHmFzRaLIr6VvAEuB64Hng08AooCtQBPQDjgW+C3QDXpT0oqQxyQ7auXhF/g2ILmPwlV2XqcxsKkGZwilAGfBDgjrdz5pZZmQ07dgLc9cDcNr4fmmNozA/j99ffjgDu3Vg1qoK7nzho7TGE7HQN6e1SaRd2xKv2024plZ2Lwc+Dwwys6+b2b/NbEnY8LzGzDaY2dtm9hszOw0YAswGzktF4M7FI7qMoX6DmqcMLsNIGhF5EHyK9l2CxYVngZ8Bwxuc41Jsx+4a3li0EQk+fmCfdIdD19JC7rzkUPIEd09dzLy129IdUn2yO9qT3VYZ0SvYpJZNE/OyRUFjT5jZ4S25kJmtIVgFdi5j1JcxSBTkR5Jdz3ZdxlkExPrI4Rzg7BjHfQ5rir32UTl7auo4fEg3+nTOjK0qRw/vwaSJQ/nbtOX8cPIcHrnmWBT+UJ9qZsai9ZGaXS9jaI36Xrue7CZco8muc7mgvs9u3t5xwT5BzWWgz6U7ANe0SAnDJ8alt4ShoetPG8Mzs9fy7rIt/OfDdZx9cP+0xLG2ooode2rp2bGIHkkctJHLIr12vYwh8eJKdiUdB/Qws2fC3/cE7gIOIqjn/a5vnHCZaO8GNaImqKUzIudiWhhuUHMZqK7OePWjcgA+MS79JQzRunYo5PrTDuCGJ+ZwxwsLOOOgfvX7E1LJxwS33d72Y57sJlq8rcd+BhwR9ftfErTJ+Qi4DvhBguNyLiFqoyao+cquy2CvS1or6R5JZ0rypbEMsnDDdjbv2EO/LiUZOQb300cOZkiPUhaX7+DJD1anJQav1227vl2K6ViUz+Yde9iyY0+6w8kp8Sa7BwIzACQVAhcD/8/MLgJuINjM5lzGqbPoPrue7LqMNRC4GRgMPAGUhy0dL4tnjLBLrreXbgKCGtl01cQ2pTA/j699fDQAf5y6GEtDx5lFG8J6XU92W00Sw8O6XR8ukVjxJrudgMhWz6OBjsAz4e/fJ+jEkDKSnpNkkm5L5fu67FM/LliKKmPwZNdlFjNbZ2Z3m9mZQG/gGqCWoJ95uaQpkq6T1KpJBpIGS3pMUoWkbZIel9TsfVvSkeFq83xJOyWtkPSgpOExzs2T9H1JyyRVSZol6aLWxJtp3l4SjOY9ZkSPNEfSuPMPHUC/LiUs2rCdqWHJRSrV99j1zWltMiIcs+xjgxMr3mR3NXBI+OszgTlmtiH8fXcgZR2QJV0WFYtzTarvs5snCnxl12WBsL3jQ2Z2GUHiez5BO7IbgZWS3pH0/XivJ6mUYADQWOBKYBIwGnhFUsdmXn4pMB74LcG9/3vA4cAMSYMbnHsrcBPBfo4zgenAo5LOijfWTGRm9Su7xwzvmeZoGleYn8eVxw0D4L43lqb0vc3MyxgSZKRvUkuKeJPdfwE/kfQYQXuxf0Q9dziwMNGBxSKpO3An3uLMxakuqvVYntfsuixjZtVm9pyZXWdmA4HjCRLXSS24zNXACOACM3vSzCYT9EMfSrCC3JSfm9nxZvYHM5tqZv8EziBY5Lg6cpKkPsC3gJ+Z2e1m9oqZXQO8QrDnI2stLt/Bxu176NWpuL41VKa6/OghdCjM5/WFG1mwrjJl77tuWxUVu6rpVlpI787pmSyXK0Z4+7GkiDfZvYlgZGUxwY3rzqjnDiGY2Z4KPydYVf5Xit7PZbnauhg1u17G4LKUmU03s++Z2bgWvOw8YLqZLYq6zlLgTYJV46beb7/Pw81sOVBOUGcccTrBZM1/NDj9H8CEWGUP2WLvqm5m1utG61payMVHDALgX++sSNn7RgZaHNivS8b/N8p0kWTXB0skVlzJrpnVmtmPzexcM7vFzGqinrvAzO5s6vWJIOkE4LPAl+M8v9GHaz8iK7veZ7f98nsB44E5MY6XAS1JmgGQdCDQB5jX4D12EwzHaPgeNPY+2fBnkw31utEuPTqoLnnig9VUVaemI+jcNWGy29/3UrZVpGZ3xead1LSjcZ/JvhfEu7IbCWaUpMslfTv8OjIhUTT/vkXAn4DbzWxBKt7T5YZ9N6gFx+o82XUZQtJISa9IWiLpDkklUc+9k6C36QFsiXF8M0E5QtwkFQB3E6zs3tfgPbba/m0ANkc9n3WypV432vgBXRk/oAsVu6p5vmxdSt5z3tqgZGLcAE9226pDUT4Du3WgutZYuWVXusPJGXElu5JKJP2F4Cf5fxCUE/wDmC/pz5KSXaTzHaAD8ON4X2BmjT5c+1HfZzcPL2NopzL8XvB74HHgUwSb0V6UFNnhU5i2qBp3F3AccIWZxUqgWyTD/2xYvmkn67ftpntpYVZtvLrkqGB195EZK1PyfvVlDP29E0Mi1Nftbmg/pQzJvhfEu7J7O/AZ4EfAKKBz+PUmgo0Sv0xINDGE7XFuAP4PKJbUTVK38OnI731OvIsp1ga1uvbzyZDLfH3N7Hdm9p6ZTQJeAF6Q1BlIVMa3hdgruI2t+MYk6WfAl4DPm9mUGO/RTft/5hhZ0d1MForur5uXhqlkrXX+IQMpLsjjzUWbWL01uauDO/fUsHTTDgryxKgs+oEgk43o5b12Ey3eZPdS4GYz+4mZLTGzHeHXHwO3kNyhEiOAEoKV5C1RDwh2/24BJiTx/V0Wi96gVuAruy7zdIj+jZndDDwLTCHob54IZQQ1tQ2NA+bGcwFJNwDfBb5mZg808h7FQMPStkitblzvk2nq63WzpIQhomtpIace2BeAp2etSep7zV9XiRmM6tOJ4gJfd0qEyMhlbz+WOPEmu8VAY/VjbxPswk2WmcDJMR4QJMAns/+mCOcAqA3z2rwGE9Qy5WNS1+4tlHRK9AEzuw14juDTs0R4CpgoaUTkgKRhBG3MnmruxZK+BtwG3GBmdzVy2nNANcEngNGuIOigk9rGrwny9tLs2pwW7bxDg/kjT81MbrK7t4TB63UTJbJJzZPdxCmI87wXgdPCrw2dRtD3MSnMbCvwasPj4adly81sv+eci6iL2qAW7OwEM6gzyM+eTyVd7ppEjHIFM7tZUqJaOt4LfAWYLOnG8P1uBVYSbPwFQNJQguEVt5jZLeGxS4FfEySzL0uaGHXdbWY2N4x3g6Q7gO9LqiSYrHkJcApB67Oss2rLTlZv3UWXkgLG9su+RO6kMb3pXFLA3LXbWLi+MmmTzSLJ7jhPdhPGe+0mXrzJ7h3AA+G0nUeB9UBf4NPAWcAV0asGZrYk0YE61xrRZQwQJL01ZtTWWf0x59Il/GG+MZWSjiMo42r4urgXGMxsR7h6fCfwACDgJeAbZhb9r6mAfPb9xO+M8PgZ4SPaVOCkqN/fAGwHvg70AxYAnzazZ8hCkRKGo4f3yMp7RXFBPmce1I9HZqziqVlr+OZpY5LyPh+uqgBgvHdiSJh+XUooLcpn0449VOyspmtpJu5VzS7xJrtTw6/XAddGHVeD5yOSXrhjZtl393EpV99nN9w3k5cnqLP6485lmnDh4EHg6Mih8KuFvzZaeI81sxXARc2csyzqvSLHrgKuivM9agnKHW5rSWyZKttajsVy/qED65Pd6z9xQML7F1dV1zJ37TYkOHhwt+Zf4OKSlyeG9+pI2ZptLN64ncOHtKhDoIsh3mT38yRuZ7BzKRNrZTf6uHMZ6M/AEOAbwHxgT3rDaZ+yuV43YuKInvTuXMzyTTuZtaqCQxOckJat2UZ1rTGmb2c6FcebTrh4jOjdibI121hSvsOT3QRo9G+npCIz2wNgZvenLCLnEijSeSE//GDWe+26LHAUcJWZ/TvdgbRX6yqqWL5pJ52KC7K6FjU/T5xzcH/++uYyJs9cnfBk94MVQWOkw4b4qm6iRdqPed1uYjTVjWGjpEckXSYpe/9vd+1aZINapIyhPtmt9WTXZaxV+GpuWkVKGI4c1p2C/BYNGs045x86EICnZ61N+PjZD1YGJeeJTqJddPsxT3YToan/i78E1AJ/BMolTZF0naQBqQnNubartUaSXV/ZdZnrJ8B3ww3BLg2mR21Oy3aHDOrK8F4d2bh9N28s2pjQa89cESS7h/nH7AlXP1jC248lRKPJrpk9ZGaXEYywPJ+gJc2NwEpJ70j6gaRxjb3euUwQmZYWSXL3TlHzZNdlpnBow1RgmaSnJf29weNv6Y4x1+XC5rQISVwQru4++cHqhF13bcUuVm/dRafiAp+clgSR9mPLNu1I+Ip8e9Ts5zNmVm1mz5nZdWY2kKAR+UsE/SHnSFog6ReSjk12sM61VH03hsgGtfBvvK/sukwl6Srg+0A34HDgxBgPlyQbKqtYUr6DDoX5HDyoa7rDSYgLDwuS3efL1rNjd01CrvnmosgPBNnZmi3TlRYV0L9rCdW1xqotyR353B60uBjJzKab2ffN7ECCUZB/Ibj5vp7o4Jxrq9qooRLRX70bg8tgNwNPAL3NbKCZDW/wGNHcBVzrvRN2YThiaHcKs7xeN2JIz1KOHNqdXdW1PF+2LiHXfCssiThuVK+EXM/tb2TvYMV80Qav222rRv9PlvRZSaVNvdjM5pvZz83sWGBQwqNzro32ruwSfo2UMaQrIuea1RP4QzMDJ1ySvLU4WLE8dmT2lzBEu/DwYHX38ffbXspgZvX/nY4flVv/nTLJ6L5BsvvRhso0R5L9mvqx9X5gnaT7JZ3c3IXMLDE/LjqXQPut7PoGNZf53gAOTHcQ7VX9imWOJbvnTBhAcUEebyza2OYd/vPWVrJuWxW9OhVxQJ/kjCF2MCYc8bxwva/stlVTye4nCD5K+yTwoqTlkm6TdEBqQnOu7fYbKhFJdn1p12WurwNXS/qMpJ6S8ho+0h1grlq9dRfLNu2kc3EBEwbmRr1uRNfSwvqNan+ftrxN13ouLIX4xLh+9Z+WucQbHSa7C9b5ym5bNdWN4SUzu5JgxvlVBHPOvwfMkzRN0rWSvLmey2iR0ty8/SaopSsi55o1D5gA/B3YAFQ3eHgP3iSJrOoeM6Jn1vfXjeXK44YB8Nh7q9jeho1qz88Jkt0zDuqXiLBcIyJlDIvLt/s+kzZqdr6fme0EHgAekNSfoAvDFcAfgDslPQP8zcyeSWqkzrVCo2UMfuNwmesWfDx7WkTqUHOthCFi3IAuHD28B+8s3cxf31jKVz8+usXXmLtmGwvWV9KlpIBjR+Tmf6dM0aWkkAFdS1hTUcWKzTsZ3stbb7dWi4ZZm9la4BfALyQdClwLXA1c2NJrOZcKe8cFN+iz6zW7LkOZ2U3pjqE9CjZdBSu7x+dwh4FvnDqay+99m3teW8JnJg6lR8eiFr3+genLAPjk4YMoKsi91e9MM7pvZ9ZUVLFgXaUnu23Qqr+pkk4hqCu7DBDBwAnnMk6j44J9ZddlEEnHpTuG9m7+ukrWb9tNr07FHNA3d4ckHDeyFyeO7kXl7hp+9FRZi167cftunggHU0w6dmgywnMNjOkX2aTmdbttEXeyK+lAST+VtAJ4AbgA+BdwvJmNSVaAzrXFfiu73o3BZabXJa2VdI+kMyW1bLnNtdkLc9cDcOqBfZBye9PVjy+YQGlRPk/PWsOfX18S9+tuf34BVdV1fHxsn/oesC65RofT6RZ4stsmTSa7knpL+rqkGcAc4FvAhwQruv3M7Fozm5aCOJ1rMTMjktNGNgznh199XLDLMAMJhkkMJuiCUy7pEUmXSeqS3tDahxfnRZLdvmmOJPmG9CzlxxceBMBtz87ju4/NZsG6yiY/8Xpkxkoeenclhfnie2eOTVWo7d7elV1vP9YWjdbZhhvPTgvPmQN8B3jQ++m6bFFbX8JA/UpNQThdosaTXZdBwvvq3cDdkjoDZwPnA38EOkiaSpAETzazNemLNDetq6hi9qoKOhTmc8Lo3K3XjXbhYYOoqq7jR5PLeHjGSh6esZKi/Dx6dCyiS4cCunYopGuHQrqUFLK2ooppS4LNezeePa6+JZZLvlHhyu6Sjduprq3Lmal+qdbUprKjCTou/M3MPkhRPM4lTMMSBtg7Sc1Xdl2mMrNK4CHgIUmFwMcJEt8bgbskvQc8YWY/TWOYOWXK3GAN58TRvSgpzE9zNKlz2dFDOHJod+59fQmvfbSRdduqwsf+55YU5vGDsw7ks8cOS3mc7VlpUQGDe3Rg5eZdLNu4w3/QaKWmkt0BZtb6RnzOpVlkbkReVP2dT1BzmUjSODOb2/C4mVUDz4WP6yRNJNgvMQnwZDdBIiN0zz64f5ojSb3RfTvzi4sPAWDnnhq27KymYmc126qqqdhVzbZd1XQoyueEUb3oVuql5Okwpm9nVm7exYL1lZ7stlKjyW7DRFfSQOCbwP8APYDzzGyOpG8A08zs7aRG6lwLxVzZlXdjcBlpjqSNBKOCpwKvATPN9v2pzMymA9MJBvy4BFhcvp2ZK7fSqbiA08a17yEJpUUFlBYVMLBbh3SH4qKM7tuZF+dt4COv2221uHrjShoPvA7UAtOAw4DIj3hDCUoeLk9GgM61VsOBErA38fU+uy7DfBU4MXxcQDBUYpukNwkS39eAd82sNn0h5qaH3lkBwFkT+tGhqP2UMLjsMSZczZ2/NkZ9iYtLvIMgfkUwwvJ0oIp9x1W+Bfw8wXE512b1PXajVnZ9XLDLRGb2e+D3AJJGAR8j+BTtROAsguR3p6TpwFQzuy1dseaS7btreOidlQBMmjgsvcE414hxA4KGLPNiFVO7uMSb7J4AXGZm2yU1/NF3PdC+P/txGSn2BjUvY3CZzcwWAYuA+6C+hOxjwKeBc4FTAE92E+Cfby+ncncNRw/vwYRBXdMdjnMxjejVkaKCPFZu3kXFrmq6dihMd0hZJ95kt6l1sF7ArgTE4lxCNZyeBlDgya7LEpKGEKzuRh4HANsJSslcG1XsquYPrwbDP687aWSao3GucQX5eYzt15nZqyqYv3Ybx4zome6Qsk68DdveAT7XyHOfBt5MTDjOJc7eld29x3yCmstUkg6Q9EVJf5e0FFgG3E6wIfiPwFFANzM7PY1h5ow7pixg685qjhneg5MO6J3ucJxr0rj+QSnDXK/bbZV4k91bgXMlTSFoeWPAqZL+BlwI/DhJ8QEg6WJJ/5a0XNIuSQvC0cXeg8M1KuYGtfDX3mfXZRJJawn2RXyHYCPwrcAYM+tnZheb2W/M7D0za1W1uaTBkh6TVCFpm6THw5XjeF77E0lTJG2SZJKuauS8V8PnGz6+0ZqYk+mNhRv527TlFOSJH547LufHA7vsF6nbnbvGk93WiCvZNbOpBDuEhwN/AQT8jHDncArajn2L4B+AHwBnEKxyXAe8IMnHibiY6vvs5sXos+vJrsssfQnKweYBZeFjaSIuLKkUeBkYC1xJsGAxGnhFUsc4LvFVoAPwTBznzgaObfB4qBVhJ03Frmq+/dgsAL5x6mjGD/BaXZf5Iiu7ZZ7stkq8NbuY2bPAs+FO4T7AJjNbkLTI9nWumZVH/X6qpM3A34CTCG7kzu2jyT67XsbgMks/9tbmXkGwmFAl6W2Cto+vE/Qz39mKa18NjCBYKV4EIGk2sBC4Brijmdd3NbO68N7/2WbOrQx7AWesm58uY21FFYcO7sa1H/NaXZcdxvbvggQLN1Syp6aOogJf52uJFv/XMrNFZvZWChNdGiS6Ee+GXwemKg6XXWL32Q2+ehmDyyRmtsHMHjOzr5nZoUBP4DLgPYJPs/4LbJX0tqRftvDy5wHTI4lu+H5LCfZanB9HbDnTqO/5snU8/v5qSgrzuOPTh1CQ7wmDyw6digsY1rMj1bXGog0+XKKlGv0/XdInW3oxSf3DcZap8LHw67xGYmn04dqHyOCIfcsYgr/yvrLbfmTjvcDMKszsaTP7jplNJFjx/Q9wJHB9Cy83HpgT43gZMK5tke7nsLAuuFrSbElfaOrkVP7ZVOys5sYng/8M3ztjLCN6d0r4eziXTLm8SS3Z94Kmfqz9naSZkq6V1KOZIE+UdA9Bb8iDExJZ0+83ELgFeNHMZiT7/Vx2ampl12t2XaaSlCfpSEnXS3oyHCP8JsEKbTnw7xZesgewJcbxzUD3tkW7j9eAbxDEeTFBmcSfJd2YwPdotZ/+dx7llbs5cmh3PnvssHSH41yL+Sa11muqZnc0wcawWwgS33nALIKb7W6Cm+QIgpWGrgQ3uk+Y2VvJDFhSJ2AyUEPj7dAwX7lr92qbnKDmfz/ai6buBZmyuivpBPbW7B4LdCLYCLwKeA6YCryWyvKxljKzHzY4NFnSE8ANkn5tZvt99pqq+/T8ddt4eMZKCvPFzy6asM89wblsEUl2P1y9Nc2RJF6y79ONJrvhRohbJP2MoL3Y6cBEYABQAmwC5gO/AR42s/ltjqYZkjoATxMk2R8zs1XJfk+Xveqa6rPrya7LLK+FXxcDj4W/n2pmyxJw7S3EXsFtbMU3kf5F0MlnAmkchnH78x9hBp85Ziij+njHSpedDh4YdA6Zs3obNbV1XnPeAs12YzCzPcDD4SNtJBUS/CNwJMEK8ofpjMdlvib77PrKv8sslxOs3K5JwrXLCOp2GxoHzE3C+8WStv/hZq/ayovz1tOhMJ8vnzwqXWE412Y9OxUzpEcpKzbvZOGG7RwY1vC65mXFjwVhL90HCWbCX5DprW1cZogktFKsPrtpCcm5mMzsoSQlugBPARMljYgckDQMOD58Lpk+Q9A/OG2LE395I2hXPOnYofTuXJyuMJxLiEMGdwNg1srcK2VIpqxIdoHfA58CfgXskDQx6jEozbG5DBWpVNinz26er+y67BFuVot+tKZ47V6C0cOTJZ0v6TyCfQ8rgT9FvddQSTWS9qm9lfQxSRcTtEADODKcanlx1DknSnpW0hckUcXWJQAAIABJREFUfVzSJyVNJtisdrOZ7WhF3G22obKKZz9cS55g0sSh6QjBuYQ6ZFBQyjBrVeKS3arqWqpzfAUo7qESaXZm+PWG8BHtZuCmlEbjskKsMoYCr9l1GUhSP+A+gv0Pfw+P5QN7Gpy6XdIBZrY+3mub2Q5JpwB3Ag8QbHx7CfhGg01jAvLZfxHkZva2egT4cviIvAZgbfi6W4BeQDXBNLXLzexf8caaaA+9s5LqWuO0cX0Z3KM0XWE4lzCHhiu7M1dWtPlatXXGrc/M5YHpy+lQmM93zxjDpBztVJIVya6ZDUt3DC771NV3Y9h7LDJBrcaTXZdZ/hc4nKBlVzQRrMyuCX99CXAtQQIaNzNbAVzUzDnL2Ju8Rh8/KY7rL2LvokRGMDOe/GA1AFf4qq7LEeMHdCU/TyxYt42de2ooLWp9GvebFz/i/reWAbB9dw3/N7mMjsUFfPLw3PvAPFvKGJxrsVjjgiO/9glqLsOcAdxrZrsaHDfgT2Z2s5ndBNwFnJXq4LLR/HWVLNm4gx4dizhuZM90h+NcQnQoymdsv87UWdCVobVWb93FH6cuRoIHvnA0N50bzJf50VNlbKisSlS4GaOpCWpNDpJocO6liQnHucSp77Mba4Oa1+y6zDIGiNWjvOFK60fhua4Z//lwLQCnj+/nLZpcTknEJrW7X11Mda1xzsEDOHF0b648bhgnjelNZVUNf3hlcaJCzRhN3QFekNRsXwtJVxHUgTmXUepirOxGEl9f2XUZpgTYZ+iCmdUC/QmG+URUhee6JpgZz4bJ7lkT+qU5GucS69BBkbrd1iW7VdW1PPZeMKbgq6cE7fgk8b0zxyLBP99eweqtDT9kym5NJbvDgefCiWUxSfoSwaaK5xMdmHNtFdlc6uOCXRbYQDAsZx9mtj5MeiOGE0yxdE34aP12lpTvoHtpIceO8BIGl1sOHxokuzOWb27VFMI3Fm5kV3UtEwZ25YC+e4esjO3XhXMOHsCe2jrumZpbq7tNJbunEzQdfzacXLYPSV8D7ibo03hhcsJzrvVijQuOrOx6GYPLMG8Ak+I477PAm0mOJes96yUMLoeN7N2JHh2LWL9tN8s37Wzx61+YGzRzOW1c3/2e+0o4eOWhd1dSXrm7bYFmkEbvAmb2LsFGiMOApyT9f/bOOzyO6urD71Evli3Lvcq9gm2MMQYM2DQTegslhJqQhDQgISEEvsQhECAhhJAKCcGhBBLA9GLAYBtww713y0UusmVLVm97vj9mZrVarZpVVrs+7/PMo907987c2RndOXPmd8/xR+MWkZ8ATwCvAFepakVrd9QwmopfxmChx4z2z5PAWSLymIjUml4tInEi8jgwBSdFu1EP7/klDL3C3BPDaHlEhIkDnGlVi7cfalLbKp/y8XrX2B1dW+IzvGca547qQVmlj2fchCzRQL2PvKo6H7gYOBWYKSLxIvJL4FGcjGbXBb1iM4x2gz/ObohoDGbsGu0JVV0A/BS4C9gtIs+LyEPu8jywG/ghcK9b16iDTfsL2JJTSHpKPKdYFAYjSpk40DF2FzXR2F228zC5ReVkdklhWI/QKlUvrfYLC3eQXxwdvswGA7Sp6lwRuQxHrrAWGIyj0/2WHo1YxDDaCM+zG2MZ1IwIQFV/LyLLgHtwYuJ6E9FKgXnAb1X1k3D1L1LwvLrnjepBvEkYjCil2tjNbVK7D9fuA5z/j7oSMo7rl87kIV35fMtBZszP4o5zhjavs+2A+kKPDfIWYCvOADwEeBd4BBgYVMcw2hXVGdSqyzxJg3l2jfaIqn6qqucDaUBPd0lT1fPN0G0cJmEwjgVG9upIWlIcuw+XNDpygqry4bq6JQyBeN7dZ+dvp6issnmdbQfU99i7BdgcsPzBLb8IJ9bj5qDFMNoVISeo+WUMYemSYTQKVa1S1Rx3MalYI9mSU8Cm/YV0So7ntCFdw90dw2g1YmOEk/y63cZ5dzftL2RHbjFdUhMY379zvXUnDcpgfP908oor+OPsyDfx6pMx3NJmvTCMViDUBDXvs8kYjPaEiLwF/FJVlzeyfhJOiuFiVf17q3Yugnh3VfUrWpMwGNHOKYO68MmGHD7bdJDLT2g4xa8nYThnZI8ac1lCISL84uLRXPHXL/jnZ9uYMrwbpw6O3AfIOo1dVf13W3bEMFoaf5xdm6BmtH+ygIUisgJn8u/nwCpV9b8/FJHewEScScNXAHswp0QN/BKGMSZhMKKfM4d346H31jNv8wF8Pq3xFjMU1RKG2iHHQjGuXzrfPnMwf5uzlW89t5THvjqGaaN71qn1bc80OEHNMCKVqhAT1MzYNdojqvpDEfkjcCcwHegEqIgcAcqAdCABJ33wYrfeCyZxqGZLTiEb9xfQMSmO0yLYA2UYjWVo9w707pTEnvxS1u45wvF9O9VZd09eCauz80lJiG2SxOcn5w0n+3AJb63cw3deWEaf9GROHpTB+P6dmTykKwO6prbEobQ6ZuwaUYvPF0LGYMau0U5R1a3AD0Tkx8ApwMlAb5yoDLnABmCequ4IXy/bL++7Xt1zR/UkIc4kDEb0IyKcObw7Ly3eydxNOfUau14iiTOHdSMpPrbR+4iJEZ64ZhwnZnbmL59uITuvhJnLspm5LJsYgZtOHcDPLxjZ7mVDZuwaUUuoOLuWQc1o76hqOTDXXYxG4mVNu3BM/bPMDSOaOHNYN15avJNPNuTw/bPqDhH24TpXz95ICUMgMTHCTacO4IZJmazdc4RlOw+zOOsQH6zZx7NfZLH/SCl/+dr4di1vaN+muGE0A3+c3RCeXZ95dg0jatiSU8iGfQWkJcVZFAbjmOL0oV1Jio9h2c68OkOQHS4qZ+G2Q8TGCFOHdz/qfcXECMf37cRNpw7gL18bz6vfOYW0pDjeW72v3WdbM2PXiFqqPbvVZd5n8+waRvTwxvJsAL5yXE8S4xr/itYwIp3UxDjOHul4a99euSdknVlr91HlU04b0pX0lIQW2/cJ/Tvz2FfHAvC7WRsbHe83HJixa0QF+SUV5BaW1SgLNUEtpo6kEu+s2sM/5m1r5V4ahtHS+HzK666x25jwS4YRbVwytjdQt7HrSXwuaoVEK9NG9+SiMb0oq/Tx6PsbWnz7LYUZu0bEo6pc+bf5nP34XA4UVBu89U1QC4yzu27PEb7/n+U89N56th4obKNeG4bREizZcZjsvBJ6dUriZDeFqmEcS5w5rBtpiXGs3XOEDfuO1Fi3J6+E+VtziYuRo9LrNoZ7LxhJQmwMb6/aw7Z2eg81Y9eIePbml7Ilp5C84gqeX5DlL68vzm5lVbWx+9nmA/7P6/fWHCgMw2jfeF7dS8f1aTDOqGFEI0nxsVwxvg8A//yspnb2xUU7qPIp5x/Xs0UlDIH0SU/mivF9UIV/fNY+35CasWtEPGv3VBuoWbnF/s9VoSaohciglltU7v9sxq5hRA6lFVW8u8p5devd7A3jWOTWyQOJEXhzRTa7Dzv3wcKySl5avAtwQoS1JredMQgReG1pNjkFpa26r6PBjF0j4tkQYKAeCjBcfSFCj4WKs3swQPqwJad9voIxjh1EZJKITBeRD0RklYhsFpEFIjJDRG4RkfqT2h9DzNmYw5HSSkb16siwHmnh7o5hhI3MLqlcPLY3FVXK/W+socqnPP7hJg4VlTOuXzoTMlt32BjcrQPnjuxBeZWPFxfubNV9HQ1m7BoRT+AM0IMBk9Q8z26NOLuesRswP+1ggIG8N7/pT6RzNubwypJdTW5nGIGIyE0ishqYD9wFpACbgUXAYZwkE/8Esl3Dd2DYOttOmLnMkTCYV9cw4N6vjKRjUhxzNh5g8qOf8K8vthMbI/zy4lFtEgP3ltOcIenFRTspq2xfyR3N2DXaBarKzc8u5ozffsqq3Xn+8ooqn/+VTF0EGru5ITy7IWUMAZ7dwCgOe+oJnVJe6WP6W2v5/YcbKa90BMGLtx/i5me/5CevrmLB1tx6+2kYdSEiq4BHgPeAE4F0VT1DVa9U1a+r6gWqOhLIAG4DugPrROSa8PU6vOQVl/PpxhxiBC52Z6MbxrFMz05J/OPGCaSnxLM3v5SE2BgevuJ4TujfNi+DJg3KYETPNA4WlvGeGwGivWDGrtEuWLjtEHM2HmDnoWJ+N2sj4OjxrvjrfCY/+qk/1WEoAr2xh4vK/YZs6Di7tWUMuYXVBvLBwvI6n0g/WrefGfOz+NMnW3h/jfOP/NbKbP/6lxa33KsbVWVLTiEl5e3r6dhoNZ4BBqrqPaq6XDV0IGhVzVfVF1X1AmASkBeq3rHAu6v3UlHlxA7t0TEp3N0xjHbByYO6MO+nU3nhGycz76dTuXpCvzbbt4hws6sNnvFFVpvttzFEjLErIv1E5FURyReRIyIyU0T6h7tfxyJLsg5x3dML+e+XLWfczdmY4/88f2suecXlLNp+iNXZ+YAjug+Fqtbwxlb6lCOlFUDoCWoxISaoHSp2jN2MVGem6r46pAyLtld7bmevz/H31WPFrpazO/752XbOeXwu5/5hbg0dcntlX34pNzyziJufXUzOkfY3OaG9o6p/VNUm/XCqulJVZzWmbnPGTxH5jYh8KCK5IqIicnM9dW8TkQ0iUiYiG0XkO408nCbz+jIvtq5JGAwjkI5J8Uwe2pWendr+IfDScX1IT4ln5e58lu883Ob7r4uIMHZFJAX4BBgB3ATcAAwFPhWR1OZsu7LKR2FZ5VG1nbvpAD97bVWzZvAfKirn5mcXc+GTn9WKj3e0rMnO5+N1+/2v2luS0ooqvvX8UhZsy+We11azdk9+rTol5VU14t02hvX7Cvyfq3zKwm2HmL/loL9sfh0SgSMllRSXV5GaEEu/jGTASTABoSeoxcW6ocfcdZVVPsorfYjAoK7OpVRXFpiVu6uP9YstB9l/pJRtB4pIjIshPlbYeaiYwwGGaVllFW8sz2b+1oOhNlcnJeVVPPHxJgB2Hy7hqblbm9S+MZRX+tiw70iLeY7ve301n20+yJyNB/jFm2vrrbvrUDEfrdvvP08eqsoDb6/jhAc+5N6Zq6isatr1W1RWyfMLd/De6r3U4Rg9JmmB8fMHQDLwTgP7uQ14CngNOB94BfiriNx+9L0Pzc7cYpbsOExyfCzTRvds6c0bhnGUJCfEcu1JznP0jPlZ4e1MAHHh7kAjuQ0YBAxX1S3g17htBr4NPB7cYO2efEb37lTvRlfuyuPWGV+SX1LB/ReO5ObTQs/3eHNFNp9uyOGskT38mUrW7TnCLc8uxqfwwdp9zL17Kp1S4v1tCssqWZudz/CeaaSnJFBZ5SM2RmqJxJ+cvZk5G504r999cRmz7jyD+NijfwZ5fflu7vrvSgDG90/n37dOJC0pvoFWjeeN5dk1PI2vLNnN6Euqf+c9eSV89e8L2JtfwvRLRnPjKQNqtL/z5eV8sHYfd583nG+ePshf7kVUuGJ8H2Yuy2bhtly+zDrkX3+oqJz8kgo6Jdc8Fs8w7dM5GcH5bYtdAy7kBLUgzW5xhVM3NSGO3unJsOMwe/NCO9gOBHgsc4vKeX7BDgBOGdyF/JIKlu/MY8O+Ak4Z3AWAP368mb/O2YoIvHzbJE4e1KX6ePcd4dH3N9A/I4WffWUkyQnVKU4/23yAogAjdObybH4ybThxDVwXBwrKEIGuHRLrrfdl1iHueGk5e/JL6ZwSz1M3TGBiI4PxF5dXkltYTt/Oyf5refP+AmZvqPbMf7B2H1tyChnSvQPgRLh4bkEWxeVVFJdX8sGaffgUenZM4vXvnUqvTs5Dypsr9vCvL5wYkS8t3kXvTsn84OyhAGzcV8CKXYc5MbMzQ7rXnnXv8ym3zviSRduda+YHZw3hx+cNr/M4qnzK31vhIaIlEJGuwEQgE8chsR34QlVrP1k2jiaPn0F0UlWfiAwBbqyjz3HAQ8DzqnqfW/ypiPQGfi0i/1TVilBtjwbvTc95o3uQmhgptzHDODb4+qT+PD1vK++u2svPLxjZLmRGEeHZBS4BFnoDNYCqbge+AC4N1eDapxfWyORRVllVw4vl8yn3zlxNblE5lT7lV++sqzExyuMf87Zxx8sreGPFHn740nLeWJ6NqjL97bV4ss+84gr+vSDL32b7wSKmPjaHa55eyPhff8TkRz9h2P3vc/JvZtfQnh4oKKuh89x2oIh3VtVM97c3v4Q3lmezM7f+SVrgeOsem7XJ/33Zzjy++e8llFaE9t6tyc7nuqcXcuKvP+LGfy3m/dV7a6XRDURVedbV4dx4SiZQMyEDwF/nbCE7rwSfwoPvrq/hJV24LZc3VuyhtMLHox9s8MsFcgvLyCkoIzUhlqtOdNJ9frBmH+v2HiEhLsbvsQ01Uc2TMPROT/YbjMXljqfecwxKiAxqniHsXRPJCbGOsUvoSWqqygF3Itt5o5wsNH/+1LkcTxnUhcHdHMPOy8BWVlnFv92nWlX4w8fV56WwrJKb/rWYTzce4N8LdvDoBzVTLM7d5PymPzp3GAO7pnKgoIwvAjzbPp/WeBtRVlnFz15bxUkPfcyEBz/m7ldWhvTYqir//XIn1/9jEXvyS0lJiOVwcQXffXFpjSgWdfHphhwm/WY2p//2U7774jL/tfLqst0AXDexH9ee5OjD/rPIua535hZz2V++4LkFO3h16W7eW72PGBG6pCaw70gp099a6+/b02665mlulp8nZm9m6Y7D/Gn2Zi548jPueW015/1hHn+bU9tIfWf1Xr+hC/DXOVtZk11tG+YVl/sf0iqqfPzofyv82vB2SA7wNvAX4E/u5xwReUVEhh7F9po8fgaiqo1xsZ8CdANeCCp/HugCTG50bxvuT0B6YJMwGEZ7o2/nFM4b1ZNKn/LiovYRhixSjN3RwJoQ5WuBUaEaFJRWcuqP/o7ExJI68nSG/mQmI37+Nmnjzgfg7VV7WLf3CL06JXHjKZmown2vr6lh7M1ev5+H318POOn4AO6duZpH3t/A4u2HyEhN4K/Xjwfg5cU7/W3vf2M1BwrKSE+JJzZG2H3YMf5yCsr43n+W+V/9/+uL7ZRV+jhnZA8evuJ4AGbM3+Hf/5rsfM57fB53/ncF5zw+t0be6/1HSmtpS99Ynk12XglDu3dg7k+m0KNjIou2H+KCJz/jhYU7yCuu9sh+mXWIa55awIJtueQWlTNv0wFuf3EZ5zw+l/8s2snq3fnM23SA15buZp2btGHB1lw27i+gW1oiP/vKCJLiY9h6oKiGEfHGcqePw3ukUV7p45mAbC5/+mSz/3NFlfKBO8lroythGNYzjfH9O5MQF8O+I6WowsQBGQx1PXm7DtU2QvfkO2W9OiWTmugZu46hFzJdsNScoFbkGo0pCbH0Tk9yt1nbs5tfUkFFlZKWFMeU4d1rrDtlcG1jd032EYrKq+jZMYnUhFgWbjvEpv3OcT7z2Xb2H3G8sOBMbAt8pb90x2H/dr2b+RvuzX1vfgnTnpjHcb+cxXVPL2Tmst1c+/RCXv5yF7ExQnys8OrS3Vz9lONd9/mUFbvy+O0HGzjn8bnc89pqyqt83DApkxW/OI9TB3fhYGE5985c7X/1r6o88/l2zn9iHtc+vYBXl+7mxUU7uO25JRwpdX6v99fsY8b8LKp86u/bleP78vVJzkPQq0t3UVJexb2vr6KwrJLJQ7ry6JXH86tLRvPp3VN4747TSYqPYdba/azenc/i7YdYt/cIXVIT+OO1J/CNyQOp8jlpoH//0SaqfMqkQRko8OgHG3j8w2pDtayyit+6DwyPXHE8t57mtL135mqKyirJOOfbjP3VLMb/+iN6fu0R+t30e95csQdfWcMPkGFiOnAtjgF5Lo439QVgGrBERMY1cXtNHj+PgtHu3+D9eJqWkPsRkTqXuli1O59tB4vo2iGByUO6Nr/nhmG0ODefNgCA/yza0agwZEczFjSFSHn/k4ETZzKYQ0DImBqVBbkk9RlJj+sfJbHXMCTGMYQyzvsuH63bz2PuzfKuc4dx4fG9+GjdflZn5/Ovz7fzjckDmbk8m5+/vhqfwg/PHspd5wzlx6+sZOaybJ5yPVA/nTac80f3pH9GCjsPFfPZ5gOkJMTxxZZc0hLjmHP3FOJjHcOtd6dkpr+1lv8u2cU9r63iuVtP9r8G/97UwYzo2ZFH3t/Ayl15rNiVx8heafz4fyspKKuke1oiOQVl/PDl5WzcV8De/FJec71pt50+kJ9fMBKfwt/c17LfnTqYzC6pvPCNk/nmc0vYdqCI+99Yw8PvreeGUwYwqGsqv3hrDaUVPi4d15s7zxnG3I05PPPFdrYfLOLnr6+u9XteOKaXX2pww6RMUhLiGNcvnYXbDrF0x2HOHdWDlbvyKCyrZHC3VB6/ZiwXPvk5L3+5kzvOHsqWA4V8sSWXDolx3H3eMKa/vY7ZG3K4+bSBbHCN3RE9O5IUH8upg7v4pR3TjuvJFtdIDOXZ9csY0pP8OuFgGYOn0wWIcR/v/DIGt25KQhy9O9Xt2fW23a1DIpMGVb/yT0uMY3TvTv4Hj60HigBYusPxMk4Z3o24WOGFhTt5fsEO7jp3mD+d4ku3TeIPH21i0fZDzN10gEvG9ia/pIKN+wtIiI3h+D6d6JGWxOMfbeKDNfu465xibp6xmG3uPhZsy2XBNsfj27tTEk/fOIGEuBi++e8lrM7O56zH5hIXKxSUVnuBO6fEc/+Fo7jS9aA/9tWxTPvDPD5at593V+/l/NE9uf+NNbz8ZXXc4IXbqj2m3zpjEBMHZPDN55bw2KyNqCr7j5QxoEsKJ2Z2RkQY2y+dlbvyuOrv81m75widU+L547Xj6BIkr7jxlAE8PW8bT3y8yX+Orj+5P0nxsfxk2nB2Hy5m1tr99ElP5uErjueMYd14c0U2P/rfSp78ZAtJCbF8d8oQ/rNoJ7sPOw95X53Qj9KKKj5Ys5fV2fmM+dWHdJxwCao+tLKcpP7OQ2VVcT45r/6q1nluD6jqAyGKXxSRe4B3cWQHZzVhk00eP48C758ieD+HgtY3G8+re/HY3g1KewzDCA8nD3TCkG3YV8C7q/Zyxfi+Ye1PpBi7Tea5H0zj1hlLSOozEnAMVp9P+fOnW7jtuSUADOvRgSvH9yU2Rph+yWi+/fxSHnpvPX/4eJPfCLphUiZ3nTMUEeGhy46ntKKKeZsOcvWEflxzUj9EhGtO6sfvZm3kv1/u8kcCuGXyQH8eas/r938Xj+LzLQdZk32E8b/+CIDThnTxx8C75qR+PD1vG89+sZ3+GSls3F/AgC4pvH/HGfx7QRaPvL/B/+o8Nkao8in/+Gw7PTom0b1jEtsPFtE/I4WLxzi64qE90vj4R2fyzqo9vLY0m8+3HKyhU/zqiX155MoxxMYIA7sO5OuTMnl39V5mLssmp6CMzinxpCXF8enGA7y7yvHCdk9L5NbJjrZ5QmYGC7cdYknWIc4d1YPP3Qllpw/txujenZg8pCufbznIC4t2+PW3N52ayVeO78X0t9exbMdhKqt8/ol5I3s5Htw7zh7Kgq25DOyayqXjevPCQuehICfEpDdPX9s7PZlN+x2vqvcK3/PexoXKoObJGCo8YzeWXq5nd29+3cZu17REBnXrwPUn9+elxTv5xcWjiI0RBrv61K1uBrYlWc49/8TMzozpm84LC3cyc9luCkorKCyr5Ixh3Zg0qAvnjurBou2H+Hjdfi4Z25tlOw+jCmP6diIpPpb+XVKYOCCDxVmHOON3nwIwomca/7hxAu+s2sv8rQcZ2j2N7581xB9N4s3vncZd/1vhPDBUOIbwuaN6MG10TyYOzKhhIPROT+ZnF4zgvtfX8LPXVvOPedtYuTufxLgYHrr8eKp8Pl5ZspvCskpuOCWT6092PLcXj+3N2yv38OC7693zOsD/BH7DpExW7srzp3H+5cWjaxm64BjOzy/Y4df7xseK3zOcFB/LUzdMIL+4grSkOH8ykEvHOfnX7/rfCn77wUb25Zf633j89PwRxMYIqYlxPHndCdwy40sKSisZ2r0Df7hmHJldUpi1dj8lFVV85biedH3ya20SaL2lUNWDIvIg8HK4+9JSNHUiYWWVzy/1MgmDYbRfRIRbThvAPa+tZsb8LC4/oU+94219Y0FLjNORYuweJrQHoi6PBWeN6MHjV4/l7ZV7mDqiOzdMcqQKnle0Y1Icj189zm/8TBvdk19ePIpH3t9AcXkV/TKS+eFZQ7nqxL7+Hzo5IZa/Xn9irX1ddWJffv/hRt5fsw9wvH3fCDHZrUNiHA9efhy3PPslADECdwdMorlhUib/+nw7b65wBnMR+O1VY0lOiOU7Zw6mV6ck/vVFFpkZKdx17jA27D3C7S8u8xscAN85c3ANYyY+NobLT+jL5Sf0ZcWuPP4+ZyvZeSVcOq4335g8sMZFFBcbw6Xj+nDpuJo3kZ25xTzx8SaKyiu5+7zhdHAnhJw4wDklniH7+WbH2D3NfbX47TMH8fmWg35tZEpCLN+YPIiM1AS/N3z93oIanl2AE/p3ZvF955CW6Bg53oSrULrSQM2uJ2MocjW7XsSF0OmCne+BMoY+fs1ubRmDp9ftlub05aHLj+en00b4JyX2z0ghLkbIziuhuLzSL0U4aUAGA7qmcurgLszf6miWY2OEe853zvvZI3vw4LvrmbMxh4oqH0vc33LCgGpH2M8uGME1Ty2gokoZ2DWV574xke5pSdw+ZTC3Txlcq6+dUxOYcctEco6UotDg5IDrTurP/K25vLtqLyt359M5JZ5nbj6J8f6HsNoRqqZfPIo12flsP1jECf3TuW5idZ3LT+jDqt15vLtqLzecksml40IH/O/aIZFvnzmIJz525C3fPmMw3YP6Gjjp0+OyE/pQUlHFvTNX85z7duT0oV05Z2S1vGTCgAy++NlZ7D5UwvCeaf7z7mnCIwURmQAMB3YDCcA3gaZGa2/y+HkUeNvpTM3+eRfyIVqABduopeUBAAAgAElEQVRyOVhYzqCuqRzfp/4JyIZhhJdLx/Xh4fc3sGp3Pst35fnvKeEgUozdtVRrwgIZBayrq9EV4/vWcJ2LwO+vHsud5wwlPSW+VpSCW04byDUn9aOwtLLWTbc+enRM4rIT+vhTV9517rCQN2mAqcO789srx/Dx+v1cOq5Pjcwm/TJS+On5w/nNe47+8KfTRtSYJR9siA7smsr3pg7mL5863trj+nTkqxPqvpmP65fO32+obaw3RP8uKTx+TW2Z4Pj+nRFx9Km5hWUs35VHbIxwsvuaf/KQrlw6rrffeP+/i0b5vY8nDchg56FiFm3P9Wt2h/esnmUfGHWhawenzcHC2vFm9/hlDMkkxzuXc7Vn17Foaxi7fs2ur0bdlIRYOiXHkxwfS2FZJUdKK+gYcH0Eyhj8fQw4x/GxMWR2SWHrgSI+2ZBDblE5XTskkNklBYBHrhjDTc8uZv+RUv7volH+SCEDu6YyuFsqWw8UsSTrMF9meUZy9XUxvn9n3vnB6WzYd4Qpw7vXikhRF429hmNihCevPYELjutFTkEpF43p7Tfq66JLh0Q+uPN0NuwtYHTvjjUesGJjhAcuPY4HLj2uwX3fcfZQBnZNJTZGuOC4Xo3qL8B1E/vTs1MS//tyF/0zUvjB2UNrPf13TIpnVO+Wi0QSJkYCT+MYuqXA+8D1TdzGUY2fR7EP3P0EGrueVrdF9uN58S8a2zuivPKGcSySFB/LdRP787c5W3nms+2Mv96M3YZ4C3hMRAap6jYAERkAnAb8rKkb65eRUue6lIQ4UhKa/rP85vLjGdo9jb6dk7loTP037atP6sfVJ4XOavKtMwZz9sgeqGrIEEvB/Pjc4WRmpHKgsIyvTezfrLBlTaVTcjzDuqexcX8BT3+2jSqfMr5/ut9IFBEev3ocV4zvS6fkeMb1S/e3nTiwM68t2+3m0PbRJz25TiPO79kNkjFUVvnYd6QUEeeBIyWh5gS1yqp6ZAwhNLsiQu/0JLYeKGJPXgkdewYYu0Ge3VCM6NmRrQeKeM6dZOhpWMF5YPjkx2eiiv+VvMfZI3uw9cA2Zq3dx0o3McWJmTUHheE902o8DLQ0sTHChQ1ct8EkxsUyNuCcHg0iUutNQmOZOrw7U4MmC0Ybqvq8iPwPR6P7HZxYuU0L3NzC42cdLHD7dT3wcUD513G8ul80dwdllVX+t2eXjG3atWoYRni48ZRMnvlsO++t2cum/QUM69F697H6iBR1/z+ALOBNEblURC4B3gR24QQxDztJ8bHcPmUwF7eAx2Fwtw6NMnTBMZyuPqkf35s6hM6u17Qt8WLKPjXXmXQVPDs6NkY4c1i3GoYuwKmDnXrbDzqTreqL81qXjGFvfik+dXTECXExpARHY/DH2a2+zD1D0wu64YUp88KW9e3sPAjtCAr15vfs1mPsjunreGsXe1KEzJrHJCK1DF2Ac0Y6obZmzM+irNLHCDc2s3FsIiK3BH5X1TJVfV9VLwXeAP7ZxE02avwUkUwRqRSRXwT150wRuQonUQTABBG5yi3z+lgB/B9wk4g8KCJTROQB4FbgF6ra7DSA8zYdpKC0kpG9OjZ6fDQMI7z06pTMNSf1Q9XJKxAuIsLYVdUiHM/GJpy4jS/iBFo/S1UL62trtC7nuTFRPS4aG1qfGUy/jBRGBHgqAyMcBONJHw4VlfujKICTrACqJwCmxNeMs1sZaoKahPbsprrGbmAihEAaY+wGG/RThners24gJw3ozMCu1YmsLrOJN8c6T4nIUhG5UkRig9atwEk40WiaMH4KEEvt+8KvcLKh/cn9/j33+ytB+/k7cDtwNTALuA74vqr+pSn9rYv3VjvqiIbenBmG0b64fcpgEmJjeHf1XjbvL2i4QSsQEcYugKruVNUrVbWjqqap6mWqmhXufh3rTBrYxe+VnTa6R5NeUdxx9lBEHN3qJWPrNvCS4mNJS4qj0qc14tFuznH+aYa6BqonPykOisYQcoKa1jR2k11D2dtW8D9kKM1uMCdmdqZvZ2eSW/+MFL/h3BAiwoOXHUdKgiML+NrJtSeEGccUE3H0ua8AB0TkHRH5o4g8BfwdxyPbJBozfqpqlqqKqk4PKp/iltdaQuznKVUdpqqJqjpUVf/a1L6GorLKx6cbnagd5x9n6YENI5LonR7g3f1kS8MNWoFI0ewa7ZSYGOHft0xk2c7DnNzIlLMeXzm+F5/+eApdOiTUSJcbim4dEikoreRgYZlfruF5Xz2j0ttGSZBmt750wWWVzkS1RM/YdY31zUGe3YON0OzGxcbwx2vH8erSbL59xqAmyVlOG9KVpfefS1J8jE28OcZR1RXAaSIyDSf6wpnABe7qXTja3WOK5bvyyCuuYECXFAYFvAUxDCMyuH3KYP775S7eWbWHH541xH+vbSsixrNrtF+SE2I5bUjXowrwPqBraq2oGKHo4kZkOBCg293sN3adf5rEOGf/ngEbyrMbF+TZ9TK7eG0DZQxe+8oqH7lF5YhUSyrq4sTMDB6+4ngGHMUNOTkh1gxdw4+qzlLVr6pqdyAd6Kqqmao6J8xda3Nmr3e8umeP7GH/I4YRgYTbu2vGrhEReJPUct3wY6rKFjeJxNAejoGa5HpnPQPWn0Et0LPrflZ1vLvlQZ7dTsnx9OiYSFmlz5+x7VBxOaqQkZLQptEuDMNDVY+oaovEqo1EZq/fD8DZI6I7+oZhRDOedvedVXvaXLtrd24jIgiOyJBTUEZBWSXpKfF0cb2tfs9uhWPAhkoqEfi9SrVaxhBX/a8w1PUUexnZGjM5zTCag4i8JSInNKF+koj8SESiXtKwM7eYzTmFpCXG1Ui2YhhGZBFO764Zu0ZE0MWfWMIxPDd7Xt3uHfyvNRODPbtu4oi4mJqXeWBEhlDGrjfJbpP75GnGrtEGZAELRWSRiPxQRMaLSI05FSLSW0QuE5FncBI3fANYFoa+tilfbHXCCk8e2pWEOLtlGUYkE+jd3ZFb1Gb7tZHDiAiqE0s4MgYvEkNgvM2keOdyLvU8u+4EtSBb1//dp0pZhafZrZ4g54VE89IYe8Zu13oiMRhGc1DVH+JkG1sMTAe+BEpF5JCI7BWREpzJaTNxspTdCYxR1cVh6nKb4aUjb+oEWMMw2h+905O5aGwvVOG/XzY5sMxRY8auERH4NbtFjuEZHIkBqg3Was+up9lthGc3vrqOl6ls474jQOOypxlGc1HVrar6A6AnTlzcnwPP4SSAeBy4GRioqpNU9d+qWhW2zrYhS3c4KbRNwmAY0cF1E53wmq8s3U1Fla9N9mmhx4yIoFuaF43B8+xWyxg8gj27VRpas+vPouarHY0BHBmDCGw7UER5pa9RMXYNowWJV9W5wNxwdyTc5BSUsiO3mJSE2BpJaAzDiFwmZHZmSPcObMkp5JMNOUwb3fqxs82za0QEXVI9GcPReHZrGrve90qfL0CzWy1jSE6IJTMjhUqfsvVAoWl2jbbmP/WtFJFjJpf00izHqzu+f+ejCm1oGEb7Q0S4ekJfAN5asadN9mmjhxERdE2rljHkFpZxqKic1IRYenVK8tcJjrMbKqlE4PcqVX/khsSgiS/D/brdI41KKGEYLcgwEXk4uFBEYkXkNiA8KYjCwBJXwnBiZucw98QwjJbkwjG9AZi9YT/F5ZWtvj8zdo2IIDUhlqT4GEorfKzYlQfAkB5pNQLMewZraUUVqhoyqQQEZlGr9gInxQcbux0BZ5KaeXaNNuYS4DYRud4rEJGvAeuBp4DccHWsrVni1+uasWsY0USf9GTG90+ntMLnTxrTmpixa0QEIkKvTskAzNt0AIAh3TrUqBMXG0NcjOBTJ8ZuZR0yhtBxdmumKx7pn6RWwP4jptk12g5V3QxcA/zNjaW7CngBKAO+qqqNjscbyZSUV7E2O58YgRP6m7FrGNHGRa53951VrS9lMGPXiBj6dnaM3U82Ok+BXua0QAKlDL66Jqj5Pbuh4+xCtYxhwdZcCssqSUuMIz2l4bTGhtESqOps4B7gMSAeuA4n1NhrYe1YG7JiVx6VPmVkr450SLS51IYRbXzleGdi2rxNBymtaN3gMmbsGhFD/4wUAHYdKgFqRmLw8FIGl1ZUUVlVR1KJmIDQYyHi7AJkdkklKT7Gbwz375JSQzJhGK2FiMwSkUeAPOAjoAB4Q9V9ejtGWLrDia87wfS6hhGV9OqUzOjeHSmpqGLhttZVZ5mxa0QMnrHrMbp3p1p1Aj27fs1ubCNkDPG1DeKhAQkrBnRJbWbvDaPRJALfBl4EzgUmAEtF5G8i8i0RmRjW3rURX7qRGE60+LqGEbWcPaI7AJ9saF3drhm7RsQQaNz2y0imZ0AkBg9/yuCKKr9mN1ZCG7uVVdXGbkKIsEbDA+J6Du5mxq7RNqjqFFXtDAwCrgR+DWwDLgT+DiwIY/faBJ9PWbbTMXZPsslphhG1THWN3dnrc2jNl1dm7BoRw8SAdKFnDO0Wsk51RAZfndEYPOPX0wglxMb4E00Ecs7I7v7PF4zp1YyeG0bTUdUsVX1dVX+pqpeoan+gG3BeuPvW2mzKKaCgtJI+6cn+iamGYUQfY/um0yU1gey8En+yqNbAVP9GxJAQF8Mfrx3H8p153HP+iJB1/J7dyip/BrXgaAyeYVtcXjt7WiDTRvfk9imDiRVhhBuKzDDCiarmArPD3Y/Wxi9hML2uYUQ1MTHClOHdeW3Zbmavz2FYj9bJlGieXSOiuHRcH6ZfMprkhNiQ6z3DtaS8ClUQoZbX1lMslFQ4gawT6jB2RYR7zh/B3dOGt1DvDcNoDEuznMlpJmEwjOjn7JGebnd/q+3DjF0jqvCiMRS5Xttgry5UyxhKykOHHTMMI7xUZ06zyWmGEe2cPrQrcTHC0h2HyS+uaJV92F3eiCo8w9VLPxgTIlxYtYzBqeNJHwzDCD/78kvZfbiEtMS4GpNEDcOITtKS4pkwoDM+hc+2HGiVfZixa0QVnrFbVNYIz25F/ZpdwzDaniVufN0TMjvXmlxqGEZ0MnW4I2X4dIMZu4bRIJ6MwfPahrpZemWeQWzGrmG0H5a4k9MsmYRhHDt4IcjmbsrB52v5EGR2lzeiilqe3RDxcz1jt8STMcSZjMEw2gtLd5ixaxjHGkO7d6B3pyQOFpazZk9+i2+/3Ru7IjJMRP4oIqtEpFBE9orIWyIyNtx9M9of1RPUGvbs+kOPxbf7fwPDaBYi0k9EXhWRfBE5IiIzRaR/I9smicjv3LG3REQWiMgZIepliYiGWC5rbD+LyipZt/cIsTHCuP7pTTlEwzAiGBFhiuvdnbOx5aUMkXCXPw+YCvwbuBj4Lk5g9YUicmI4O2a0P6o9u66xG2qCmltWbJpd4xhARFKAT4ARwE3ADcBQ4FMRaUxqwGeA24BfABcBe4FZIjIuRN1ZwClBy9zG9nXFrjyqfMro3h1JSbAw8IZxLOHX7W5s+dTBkTCavAz8RQPyyInIJ0AWcAdwY5j6ZbRDPEmC57UN5dmNj61pEJuMwYhybsNJPTxcVbcAiMgqYDPwbeDxuhq6b9C+Btyqqs+6ZXOBtcADwCVBTQ6q6sKj7egSSyZhGMcspw7uQkJsDCt25XGoqJyM1IQW23a7d2mp6kENSpisqvnAJqBPeHpltFeSXElCoWvIxsXWNnYT4pyywlLP2G33/waG0RwuARZ6hi6Aqm4HvgAubUTbCuC/AW0rcZwQ00QksSU76kVimGDxdQ3jmCM1MY6TB2WgCvM2tayUISLv8iKSARwHrA93X4z2hWe4FpQ6gakTQkxQ8zy7BZ6xa5pdI7oZDawJUb4WGNWItttVtThE2wRgSFD5xSJSLCJlIrKwKXrdKp+yfGceABMsc5phHJOcOawbAHNaWMoQqXf5PwECPFFXBRGpczGiFy9BhGfIhkoF7Dd2TcZwTGBjARnA4RDlh4CGrMr62nrrPd4GfgBMA64HSoHXReTrdW088Fwk9xpMYVklFXn76NkpuYFuGYYRTXjjwO2XOXNfX5u/HomJbbFxus2NXRE5p44Zu8HLnDra34ujIft+4Gs5w4DaMob4EJ5dzwAuLHO8vyZjMIzmo6o/UNXnVPUzVX0VOBtYAjzcmPZJfUcDULZ7Xet10jCMdk3loWwqDu8lNqUTCb2Gtth2w3GXnw+MbMRSa+KZiHwH+A1wv6r+q76dqGqdixG9eF5av2c3lLHrlplm99jAxgIOE9qDW5fXtrFtodrDWwtVrQJeAfqKSK866viXr/9oOgBP/t8dx9K5MQyDmmPBNy+YBMDDM95qsbGgze/yqlqsqhsasewMbCciNwB/BX6vqg+1db+NyCBYsxsfFyoag7h1PM2uyRiMqGYtjvY2mFFAQ27UtcBAN3xZcNtyoLFv1+q9Y6kqi7c7dvPJg2xymmEcy0wZ7sXbbTndbkS4tETkcuBZ4J+qene4+2O0X7ykEhVVzr01pGfXNYgrfXXXMYwo4i1gkogM8gpEZABwmruuPt4G4oGvBrSNA64BPlTVsroaBtTbqar76tvJ1gOFHCwsp2uHRAZ1bUzoX8MwopVJg7qQFB/Dqt357MkraZFttvs4u26mnpeAlcAMEZkUsLpMVZeHp2dGeyRYkhBKsxtclpRgnl0jqvkH8H3gTRG5H8fL+mtgF/CUV0lEMoGtwAOq+gCAqi4Xkf8CT4hIPLAduB0YiDMJzWt7HU4Ys/fc7fYAvgeMB65rqIMLt1V7dY+hiYOGYYQgOSGWqcO78/6afby/pt7n5EbT7o1d4CwgEWfQ/CJo3Q5gQFt3yGi/BEdWiK8nGoNHkml2jShGVYtE5CzgD8DzOJFsZgN3qmphQFUBYqn9xu8W4CHgQSAdx/FwvqouC6izHegO/A5Hz1uEMzntfFWd1VAfF7kShkkDTcJgGAZccHwv3l+zj/dW722R7bV7Y1dVpwPTw9wNI0JICoqZm1jPBLXqNubZNaIbdw7ElQ3UycIxeIPLS4AfuUtdbRfiOCaajM+nLNh6EHBeXxqGYZw1ojuJcTEs3dHQHNrGYS4tI6qo5dmtR7PrkWzGrmGEjdXZ+RwsLKdPejJDuncId3cMw2gHpCbGMdWdqNYSmLFrRBXBnt3Q0RjMs2sY7YVPNjgzrqeO6GZ6XcMw/FwwJmTEwqPCjF0jqmiMZ9cLPeaRnGD/BoYRLmZv2A/Qol4cwzAin7NdKUNLYHd5I6pIDPLshkoXHFxm6YINIzxsySlgTfYR0hLjOG1I13B3xzCMdkRqYhznje7ZItsyY9eIKoKfAuvLoOZhMgbDCA+vLs0G4MIxvez/0DCMWlw5vk+LbMeMXSOqEJEantvGxNlNtji7hhEWXl26G4ArT+wb5p4YhtEeOX1otxbZjhm7RtSR2ICxGyxjsDi7hhEeDhaWMaZvJyZkdg53VwzDaIfExrTMpFW7yxtRR+Dr0FCaXfPsGkb74Z7zR1gUBsMwWhUzdo2oI9CzmxBb+yaaEBSOLMkmqBlGWLj7vGE2Mc0wjFbHjF0j6mhIxhBYlhAbQ0wLvSYxDKNpfP+soeHugmEYxwBm7BpRR0MyhsCy4CQUhmEYhmFEF3anN6KOQM9uXAjPbsekeP9n0+sahmEYRnRjxq4RdQQasB2T4mqt75yS4P+ckZrYJn0yDMMwDCM8mLFrRB2BBmxGakKt9YHGcLLJGAzDMAwjqrE7vRF1dEkN9NzWNnYDaakYfoZhGIZhtE/M2DWijq4dGm/sWnxPwzAMw4huzNg1oo6UhLiQn0Nhjl3DMAzDiG7M2DWijsYYsP0zUgA4bbAFtDcMwzCMaKZ+t5dhRCCDunVosM5/vz2JTzbkcNWJfdugR4ZhGIZhhAtR1XD3oUUREQWItuMyGo+q8sKinYzrm87xfTuFuztGO8bTbKuqCVraEBunDcNoLC0xTpuxaxjGMYsZu+HBxmnDMBpLS4zTptk1DMMwDMMwohYzdo1jAhGxMGOGYbQa0TrGRONxReMxQfQeV0tgxq5hGIZhGIYRtZixaxiGYRiGYUQtEWfsisi1IqIisjvcfTEMw4gERKSfiLwqIvkickREZopI/0a2TRKR34nIXhEpEZEFInJGiHoxInKviGSJSKmIrBSRK1v+aAzDMJpGRBm7IpIOPAHsC3dfDMMwIgERSQE+AUYANwE3AEOBT0UktRGbeAa4DfgFcBGwF5glIuOC6v0amA78GfgKsBB4RUQuaIHDMAzDOGoiKvSYiDwNZOIMtueoaq2MABbSxghFQOiSMPfEaE8cC6HHROQO4HFguKpuccsGApuBn6rq4/W0HQusAG5V1WfdsjhgLbBRVS9xy7oDu4BHVPWXAe1nA91UdUzQdqNunI7WMSYajysajwmOieOK/tBjInIa8HXge+Hui2EYRgRxCbDQM3QBVHU78AVwaSPaVgD/DWhbCbwMTBORRLd4GpAAvBDU/gXgeNe4NgzDCAsRkS5YROKBp4HfqeqWxoTWsPAbRijsujCOQUYDb4YoXwt8tRFtt6tqcYi2CcAQ9/NooAzYEqIewChge/DGo/H/MRqPCaLzuKLxmCB6j6s5RIpn9x4gEXg43B0xDMOIMDKAwyHKDwGdm9HWW+/9zdPa70+D6xmGYbQ5be7ZFZFzgI8aUXWuqk4RkSHAfcDlqlraUKNo1t4ZhmFEAzZOG4bRloRDxjAfGNmIet5rsydxZhIvdKMxgPP6TNzvZapa0vLdNAzDiAoOE9qDW5fXNrhtZh1todpzexhIFxEJ8u4G1zMMw2hz2tzYdbVfG5rQZBTOYBtqUD4M/BG4swW6ZhiGEY14mtpgRgHrGtH2chFJCdLtjgLKqdborsWRmg2mpm53lPu3of0YhmG0GpGg2b0WmBq0zAIOup//HL6uGYZhtHveAiaJyCCvQEQGAKe56+rjbSCegIlsbuixa4APVbXMLf4AJ2rD9UHtvw6scaM/GIZhhIV2b+yq6kJVnRO44CSVKHO/b2lO5h4RmeFmZAtenmjVAzNahLbIDGVEHs28LkKNBxoiiUKk8A8gC3hTRC4VkUtwojPsAp7yKolIpohUisgvvDJVXY4TduwJEfmmiJyNE3ZsIPDLgHo5OLF87xWRH4nIFBH5G3AWcG/APo76vLQHRKSviPzJHSuK3etiQIh6ETO2iMhVIvKaiOxw+7pRRB4WkbSgep1F5J8iclBEikTkYxE5Plz9bggRmSYin4jIPhEpE5HdIvI/ERkVVC/Sr8kP3OvwwaDyiDlf7ngRaszNC6p31McUEaHHGsGvgbtxJrItxfEGvyIiF6nqe41ofwAnnmQge1u2i0ZLI9WZocpwMkMp8CBOZqgxqlrUwCaeAS4EfgJsw4nhPEtETlHVFa3Xc6M1aYHrAmAGAYagy6aW7GdboapFInIW8AfgeUCA2cCdqloYUFWAWGo7QW4BHsL5DdOBlcD5qrosqN59QCFwB9AT2AhcrarvQIudl3AzBLga5z7zGXBeHfUiaWy5G9gJ/BzYDZyAkwlvqoicqqo+EREcL/8A4Ac4EsJ7cc7dOFXdHY6ON0AGznn6K849vj/wM5z5P8er6o5IvyZF5DpgbIjySDxfAD8Evgz4Xul9aPYxqWpEL0B3nAv1V0Hls4FVjWg/A9gd7uOw5ajO/R1AFTAkoGyg+w/yowbajsUZ2G4JKIvDuUG/Fe5jsyU814VbV4EHw30c0bY097y0hwWICfj8TfdaGRBUJ6LGFpwMd8FlN7rHcJb7/VL3+9SAOp1wJh4+Ge5jaMKxDneP48fu94i9JnEmne4DrgsesyLtfAFT3P6eU0+dZh1Tu5cxNALL3HPs0haZoYzIoznXhdF6RPx5UVVfI6pF1NiiqgdCFHvetT7u30uAPar6aUC7fBxPW0ScO5dc96/nMYzka/JRHD38SyHWRcv5CqRZxxQNxm5jMvc0RHdXA1IpIptE5B4RiW3RXhqtwWhgTYjytTR83huTGcqITJpzXXjc7ur8il3d3+kt171jlpY4L5FANIwtZ7p/17t/6zt3/UWkQ5v06igQkVgRSRCRoTjSpH2AZyBG5DUpIpNxvO/fq6NKpJ6vF0WkSkRyReQ/QdrpZh1TNGh2m5u5ZwWOrmctkARcjpOpbSjOayqj/dIWmaGMyKM51wU4b4XeAfbghD38CfCJiJyrzgRZ4+ho7nmJFCJ6bBGRPsADwMequsQtzsCZ5BiMd0ydcfTa7ZFFwInu5y040owc93vEXZMikoBjtD+mqhvrqBZp5ysf+D0wFziCoxv/ObBARE5wz1ezjqndGbvSxAxrzd2fqgZHXXhPRAqBO0XkUVXd3Nx9GIYROajqDQFfPxORN3E8Cg8Ck8PTK8NofVzv2Js4r/lvCXN3WoobgI7AIJzJeB+JyGRVzQprr46enwLJOJNGowJ1or4sDyiaKyLzgMU4k9bub+4+2p2xS9MzrLVG5p6XcBJVTADM2G2/tEVmKCPyaM51UQtVLRCRd4FvNLdjxzgtel7aMRE5tohIMo7+cRBwptac3V7fufPWt0tU1ZNiLBKR93G8gz8DvkOEXZPua/37cN46JwbpvxPFySpbQASfLw9VXSYim4CT3KJmHVO7M3a16RnWWjNzT7A0wmhftEVmKCPyaM51UR82HjSP1jov7Y2IG1tEJB54FcfBc66qrg6qspbQYdZGATu1Zgi7douq5onIFqp105F2TQ7CkVsGT8gHx2t9N44EICrOl4s37jbrmKJhglprZO65HucH/rKhikZYaYvMUEbk0ZzrohYi0hG4COeVmnH0tOh5acdE1NgiIjHAizgJQC5T1YUhqr0F9BGRMwPadQQuJoLOnYj0AEYAW92iSLsmV1A7o+xUd90L7uctRMH5EpEJOKHivHG3Wcckted1RR4i8giO7ODnwDKcQeXbwCXqBjR3680GMlV1iPs9EyfI+ss4F0gizgS1m4GnVPX2NjwMo4mISCpOgPsSHLKIcs8AAA5dSURBVE2P4iQYSQPGeE967nneCjygqg8EtH8ZJ3TdT4DtwO04Rs2pWjtgvhEhNOe6EJG7cQbYT6meoOaVna2qn7Xt0UQPjT0v7R0Rucr9eDbOq/Dv4iQtOKCqc906ETO2iJPp7js4GtB3glbvVtXdrkH8OdAP55i8gP5jgLGquqsNu9woROR1HHtgFc6kp2HAXTgJTyaq6qYouiYVeEhV73e/R9T5EpEXcf5PlgF5ON7pe3HkquNV9WCzjyncwYRbYsHJ+nM/sAMnDNkq4KoQ9eYAWQHfM4A33Hal7g+7DPg+AcHDbWm/C05WnNdwBrMC93wOCKozAGcQmx5UnoyT4nSfe/4XAVPCfUy2hO+6wPESfAEcxHljlIvjNZgY7mOKhqUx56W9L+41E2qZE1AnYsYWHA1rXcc0PaBeBvAvHM1xMU7iprHh7n89x3UPTqSlPLe/G3GiGAwIqhct1+SDQWURc75wjNZVOFEZKnBSmT8N9GqpY4oKz65hGIZhGIZhhCIaNLuGYRiGYRiGERIzdg3DMAzDMIyoxYxdwzAMwzAMI2oxY9cwDMMwDMOIWszYNQzDMAzDMKIWM3YNwzAMwzCMqMWM3QhCRE4Rkf+JyB4RKReRXBH5SERuEpHYcPfPQ0SmiIiKyJSjaDtdRM4KUT5DRLJaon9N6Eu625/xIdbNEZE5bdmfhhCRt0Tkz62wXRGR5SLy05betmEYhmG0NmbsRggicidOsPsMnGDZ5wC3ApuAv+Fk54kGfomTtjKYX+Nkt2tL0nH6U8vYxcma9N227U7diMgZOHnDH27pbasTjPsB4F4RyWjp7RuGYRhGa2LGbgTgGjKPA39W1XNU9XlVnaeqb6rq94DjcVLtNXc/iXWUi4gkNHf7zUFVt6rq8nD2IRBVXaeq68LdjwB+ArytqtmttP23cDJBfbOVtm8YRjtFRHqJiE9EJovIVSLymojsEJESEdkoIg+LSFor7PdJEXnH/ayNWLLcuje734eE2Gacu256E/typ4isdtPWGhGGnbTI4B6c9HghXyO7huAq77uITBSRj0WkUESKRGS2iEwMbOPKAna70oj5IlIC/NZdlyUiL4jIrSKyASgHLnTXpYjIoyKy3ZVSbBeR+xoaAETkPBF5T0T2ikixiKwRkR8Hyi/c/N4A9wUMXtMD+psVtM1eIvKciBwUkTIRWSUiXw+q4w16k0TkRRE54spAnhSRpHr6O4DqB4h/BPTnZnd9DRlDgHTjMhF5SkQOiUieiDwhIrEicpKIfO6ej7UiMi3EPs90z1WBW2+WiBxX3+/qtusNfAX4Tx3HfoaIvOFeD7ki8hcRSQ6oFycivxaRrSJS6v6en4vIZK+OqlYBr2DGrmEci1wKHADmA3cDVcDPgfNx3izeDnzUkoagiAwGvgNMd4tOCVr2AbOCylrz7d9TQDfgplbch9FKxIW7A0b9uMbgVOANVS1tRP0xwFxgHXAzTs7snwFzRWSSqq4MqN4JeBl4DGfgKglYNxUYB/wKyAGyRCQOZ3AZhSMrWA1MAv4PR17x43q6Nggnj/WfcDyEE3AGsW5u/8AZrBYAM3AGFoDddRxnqnucnd2+7wK+DjwvIimq+nRQk+eBl4Ar3P1MBw7jyBRCsdetOxNHGvCWW761nmMEeMJtcw1wBnA/EIsjO/kdkO2WzRSRTFU96B7PhcCbwLvucYDzkPOZiIxR1V317PNcdx+f1bH+BeB/wF+BicAvgFSc68Pbz13AfcAKoCPO+QmWLMwDfiAig1R1W72/gmEY0cRlOG+OfCJysaoeCFg3V0QOAf8GpgCftNA+7wRWquoSAFVdGLhSRMqAg8HlrYWqlojIczjG/rNtsU+jBVFVW9rxAvTAMVgfbmT9V4E8ID2grCOOZ3hmQNkMd7uXhthGFlAM9Awqv8Ftc0ZQ+X043t/u7vcpbr0pdfRRcB607sMxOGMC1inwYIg2M4CsgO/fD7UP4GMc4zzW/X6zW+9XQfXeATY18FsOcNt+M8S6OcCcgO/eMf8rqN4yt3xyQNkYt+ymgLItwOygth2Bg8ATDfTzb0B2iHLv2P8e4nxVAcMCfouZ9e3DrTfY3d7Xwv1/YYsttjR9wXFgVADfqGP9L93/8e4BZR2BMuCierY70m13Q0DZf9yyupYP69leontvuKueOlnAC3Ws88a+ISHWxbnrprvfp9TTxxlBbce75aeG+1za0rTFZAzRxxnAO6qa5xWo6hEcz+SZQXUrcAydUCxU1X1BZecDO4D57qvvONfb+yEQj+PlDYkrOXhKRHbgGMYVwIM4k8C6N/roqjkDx8CbE1T+Ao63eFRQ+btB31cD/Y9ivw3xftD3DUCRqn4eVAbQD0BEhuIYki8G/a7FOJ7uMxrYZ2+cV4x18b+g7y/jSJg8acuXwAUi8pA4mry69NnePno30B/DMNonfwHWAv+qY/1692+g7O0CnDH743q2691b1geU/QZYg/PGyZMZnO+uewT4Vj3bm4Rzb6jrbVVjiQ0cU91xNThy0TJqSyTuc9etD6q7Aiig+jiMCMFkDO2fXBx5QWYj62fgvIIPZh/OK/9ADqijxQxFqG10d/tRUUebLqEKXR3XWzhG0nQcY68E59XYfUCd2tl6qO84vfWBHAr6XobjPWhpDgd9L8fxtPtR1XIRgerj9oz9Z9wlmJ0N7DMJ53jqYn8d3/u4f3+DIy35Oo4kpFBEXgV+oq7MwsWTuSRjGEZE4c4TOBW4Rl03ZQi8sSbwfnMZMEvrkNGJSB+caC0fqys5AFDVNSLSF8f5stCt6xnF76hqVj3dnYTjQV1VT53GsKGhCq4zyC+FEJFhOFKF13DnsQTU9YnISupx7BjtEzN22zmqWulOhDpXRBJVtT6jBhyjrmeI8p7UNsTqGvDqWpeLM2nr6jraZNVRPhhHA3qDqr7gFYrIxfXsvyEOAcNDlPcMWB8p5Lp/7yW096S8Ee0H1rO+B443J/A7OPphVLUCeBR4VER64oSxexxIwdEee3gPEIEGsGEYkcE3gXzgda9ARG4ENqjqYreog/s3xl2fgDP59XuhNigiHXDmGlQCtwSty8TxzgbOExmHc29Z3UBfewNHVLWhsa8hLqf2vI9YAozbQESkM/A2jqzshjoeCg4Aw5rZL6ONMWM3MngERyP6W+CO4JUiMhBIUyciw1ycV9Jpqlrgrk8DLv7/9u4tROoyjOP49wHzoi4qETqAUERBCZpgB7LQik4asiWFeMCNDhYFYVBdSBCVkGZ2gIhSQbYMitBMMAm8cO1ARV101A4UIkKIBKIFaf66eN6R2dmZ2XVPzgy/D8iws///f/7rwn+f93mf93nLNYZjOzAPOCxpwBFzldPL64mMcEScBiysc+y/DC5zuBO4KyJmSPq06v0FZM3uSLQFqwwsRjuTuYccKEyW9PwQzt8N3BER4yQdq/P9u+m7aGQ+cBz4ovbAUrqyLiJmA7WdICoB9Z4h3KOZnVqzgM/K4JbSjeZN4H6gEuxWkgWVTjQ3kM/vfuVupaPLVnLx8UxJtUHl1PJanZ2dRq69ODTAvQ40WzVY30v6tfqNUsrQT/mb9H757JmS/ql3HDnD5dmtNuNgtw1I6o2Ix4A1EXEZuVhrL1mWcCM5Yl9APlSeJTNzOyJiJTmKfpJ8YD0zzFvZSI7ed0TEi+SIfTyZuZ0LdEn6u855P5G1visi4j8y6F3W4DN+BOZExHYyE71f0v46x20gA/9NEbGcHL0vJDsTLG1SnnEy/iSzpvMj4lvgCPC7pIPNTzs5khQRDwNbSiblPTJ7eg457bhX0poml+glu2ZMIevPas2OiBfI2uoryUUoPZJ+AYiILeTv8hvy/3waWZP2Rs11riJ/d2Oy+tnMRkZEnAlMpG83mevJUq59Ne8dIzcwgixh2Fm9BqRcrxIYTgduklQvUzuFLI+qHhxPo2+mt5GDZFZ4LL0GXAHMqLNepdoEPLvVdrxArU1Iehm4lqz/XE1m6jaQq2CXkiNsSnZ3FnCIbAXzFnCYHKkO5iHT7B6OArcAa8nFBdvIAHgJ2X+x7pRTmYrqIutpe8iHSi+Zsa71CBlUbiUXTtVdxCDpCLko4uNynS1kJmGx+rcdGxJJx8mBxNlkecFXZIZ8xEnaRi5EOwNYR7Z4W0VmWj4f4PRdwP4m97aInHbbTLaHW0vf3d96yd3X1pPZ+4fKZ9f2db4d+LDBgMbMWlclsVW9PuJBMhgdDxARE8lZoK2S/opcWDAX+KD6QmUNxkYy69ulxq2/pgI/VBIPJUC+lMEFu7uB8aXmd9RFxDJyR9L5DQL3ahfi2a22E43r1M2sXURuvrGQbCem8l432Q/y4tqpvCFc/3yyl/HNknYM727NbKxFxF6yJncpGYjeSQaV48gB9nKyi810SXsi4mpyoD2pukQhIl4nA+UV9C9v2Fc5NiJ+BnZJurd8PYmckXxU0qsD3OsFZCnFPEmbGhzzB/CJpEV1vtdNg2dfKWM4SrajfDoiriETBj30n806IOm3qnPPIteDPCBpXbOfwVqLM7tmneElctpv3ihd/3FyOtOBrll76ian33uA68hd0VaRtfnvklneWZIqWcsu4Os6tbi3ldflZDBc/e8+OLHpz0X0rdc9QLYiWxkRi5vdaOnU8CWjNJNW4xIyFuqm/8/zVM2xc8gZzM1YW3Fm16xDRMStwARJ75SvuxmBzG6ZznyCLGGo7TtpZh0ocqv4tyU9d4o+vxt4BTivVUqnIuIjcte2psG6tR4Hu2ZmZtZSSrnBd8B6Satb4H4uJzvYTB5uWZiNPZcxmJmZWUspbRTvIXeSbAXnAt0OdNuTM7tmZmZm1rGc2TUzMzOzjuVg18zMzMw6loNdMzMzM+tYDnbNzMzMrGM52DUzMzOzjuVg18zMzMw61v+RaXe0Kq+87gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(10,5))\n",
    "subplot(1,2,1)\n",
    "set_fig_properties([gca()])\n",
    "plot(shc['t'], (shc['Ki']+shc['Ko'])/Ly, linewidth=2)\n",
    "xlim([-0.5, 0.5])\n",
    "gca().set_xticks([-0.5, 0, 0.5])\n",
    "ylim([-4, 10])\n",
    "gca().set_yticks(range(-4,11,2))\n",
    "ylabel('K (eV/ps)')\n",
    "xlabel('Correlation time (ps)')\n",
    "title('(a)')\n",
    "\n",
    "subplot(1,2,2)\n",
    "set_fig_properties([gca()])\n",
    "plot(shc['nu'], Gc, linewidth=2)\n",
    "xlim([0, 50])\n",
    "gca().set_xticks(range(0,51,10))\n",
    "ylim([0, 0.35])\n",
    "gca().set_yticks(linspace(0,0.35,8))\n",
    "ylabel(r'$\\kappa$($\\omega$) (GW/m$^2$/K/THz)')\n",
    "xlabel(r'$\\omega$/2$\\pi$ (THz)')\n",
    "title('(b)')\n",
    "tight_layout()\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(a)** Virial-velocity correlation function. **(b)** Spectral thermal conductance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The figure above shows the virial-velocity correlation function $K(t)$ and the spectral thermal conductance $G(\\omega)$. See [Theoretical formulations](https://gpumd.zheyongfan.org/index.php/Theoretical_formulations) for the definitions of these quantities. Using the spectral thermal conductance, one can perform quantum correction, see [[Li 2019]](https://doi.org/10.1063/1.5132543)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('../diffusive/Gc.npy', Gc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. References\n",
    "- [Fan 2019] Zheyong Fan, Haikuan Dong, Ari Harju, and Tapio Ala-Nissila, [Homogeneous nonequilibrium molecular dynamics method for heat transport and spectral decomposition with many-body potentials](https://doi.org/10.1103/PhysRevB.99.064308), Phys. Rev. B **99**, 064308 (2019).\n",
    "- [Li 2019] Zhen Li, Shiyun Xiong, Charles Sievers, Yue Hu, Zheyong Fan, Ning Wei, Hua Bao, Shunda Chen, Davide Donadio, and Tapio Ala-Nissila, [Influence of Thermostatting on Nonequilibrium Molecular Dynamics Simulations of Heat Conduction in Solids](https://doi.org/10.1063/1.5132543), J. Chem. Phys. **151**, 234105 (2019).\n",
    "- [Lindsay 2010] L. Lindsay and D.A. Broido, [Optimized Tersoff and Brenner emperical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene](https://doi.org/10.1103/PhysRevB.39.5566), Phys. Rev. B, **81**, 205441 (2010).\n",
    "- [Tersoff 1989] J. Tersoff, [Modeling solid-state chemistry: Interatomic potentials for multicomponent systems](https://doi.org/10.1103/PhysRevB.39.5566), Phys. Rev. B 39, 5566(R) (1989)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
